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ABSTRACT 


Local-Mesh,  Loeal-Order,  Adaptive  Finite  Eleaent 
Methods  with  A  Posteriori  Error  Estimators 
for  Elliptic  Partial  Differential  Equations 

Alan  Teiser 

Tale  University,  1981 


The  traditional  error  estiaates  for  the  finite  eleaent  solution  of 
elliptic  partial  differential  equations  are  a  priori,  and  little 
information  is  available  froa  then  about  the  actual  error  in  a  specific 
approximation  to  the  solution. ->^ius,  in  the  traditional  finite  eleaent 
method,  the  choice  of  discretization  £V4^ft^tojthe  user,  who  must  base 
his  decision  on  experience  with  similar  equations.  \ 


v*  In  recent  years,  locally-computable  a  posteriori  error  estimators 
have  been  developed,  which  apply  to  the  aotual  errors  coaaitted  by  the 
finite  element  aethod  for  a  given  discretization.  These  estiaators  lead 
to  algorithms  in  which  the  computer  itself  adaptively  decides  how  and 
when  to  generate  discretizations.  So  fsr,  for  two-dimensional  problems, 
the  ooaputer-generated  discretisations  have  tended  to  use  either  local 
aesh  refinement,  or  local  order  refineaent,  but  not  both. 


In  this  thesis,  we  present  a  new  class  of  loeal-aesh,  local-order. 


j i 


discretizations.  We  present  several  new  locally— computable  a  poateriori 
error  estimators  which,  under  reasonable  aasumptions.  asymptotically 
yield  upper  bounds  on  the  actual  errors  committed,  and  algorithms  in 
which  the  coaqtuter  uses  the  error  estimators  to  adaptively  produce 
sequences  of  local-mesh,  local-order  discretizations. .  We  present 
theoretical  and  numerical  results  indicating  the  expected  and  actual 
behavior  of  our  methods.  \ 
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CHAPTER  1 
Introduction 


The  traditional  error  estimates  for  the  finite  element  solution  of 
elliptic  partial  differential  equations  are  a  priori,  and  little 
information  is  available  from  them  about  the  actual  error  in  a  specific 
approximation  to  the  solution.  Thus,  in  the  traditional  finite  element 
method,  the  choice  of  discretization  is  left  to  the  user,  who  must  base 
his  decision  on  experience  with  similar  equations. 

In  recent  years,  locslly-oomputable  a  posteriori  error  estimators 
(e.g.,  Babulka  and  Rheinboldt  [8])  have  been  developed,  which  apply  to 
the  actual  errors  committed  by  the  finite  element  method  for  a  given 
discretization.  These  estimators  lead  to  algorithms  in  which  the 
computer  itself  adaptively  decides  how  and  when  to  generate  discreti¬ 
sations.  So  far,  for  two-dimensional  problems,  the  computer-generated 
discretizations  have  tended  to  use  either  local  mesh  refinement  (e.g., 
Babulka  and  Rheinboldt  [10],  Bank  and  Sherman  [13]),  or  local  order 
refinement  (Babulka,  Szabo,  and  Katz  [11]),  but  not  both. 


In  this  thesis,  we  present  a  new  class  of  local-mesh,  local-order 
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square  finite  elements  which  can  easily  accomodate  computer-chosen 
discretizations.  Ve  present  several  new  local ly-computable  a  posteriori 
error  estimators  which,  under  reasonable  assumptions,  asymptotically 
yield  upper  bounds  on  the  actual  errors  committed,  and  algorithms  in 
which  the  computer  uses  the  error  estimators  to  adaptively  produce 
sequences  of  local-mesh,  local-order  discretizations.  Ve  present 
theoretical  and  numerical  results  indicating  the  expected  and  actual 
behavior  of  our  methods. 

Ve  begin  in  Chapter  2  by  considering  bisection-type  local  mesh 
refinement  with  square  elements.  Ve  present  the  class  of  "l-irregular" 
meshes,  a  variant  of  a  class  of  meshes  suggested  by  Bank  and 
Sherman  [14],  such  that  any  bisection-type  locally-refined  mesh  can  be 
converted  into  a  1-irregular  mesh  with  no  more  than  a  constant  times  as 
many  elements,  and  the  resulting  finite  element  matrix  (with  bilinear 
basis  functions)  has  no  more  than  a  constant  number  of  nonzeroes  in  any 
row.  Conversely,  we  show  that  there  are  simple  bisection-type 
locally-refined  meshes  with  essentially  dense  finite  element  matrices. 

Ve  discuss  simple  data  structures  for  handling  1-irregular  meshes. 

In  Chapter  3,  we  extend  the  finite  element  spaces  considered  in 
Chapter  2  to  include  basis  functions  with  locally-varying  polynomial 
orders. 

In  Chapter  4,  we  develop  new  a  posteriori  error  estimators.  Ve 
present  a  new  error  estimator  which  has  two  main  advantages  over  the 
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estimator  presented  in  Babulka  and  Rhe inbold t  [8]: 

1.  Under  reasonable  assumptions,  the  new  error  estimator  is 
asymptotically  an  upper  bound  on  the  norm  of  the  true  error. 

2.  The  new  error  estimator  can  be  computed  locally  in  each 
element.  (The  error  estimator  in  [8]  must  be  computed  over 
more  complicated  looal  regions.) 

Under  suitable  assumptions,  we  show  that,  like  the  estimator  in  [8].  the 
new  error  estimator  is  no  larger  than  a  constant  times  the  norm  of  the 
actual  error.  We  also  present  several  cheaper  error  estimators,  and 
discuss  their  properties.  Under  suitable  assumptions,  Babulka  and 
Hiller  [4]  have  recently  shown  that  the  error  estimator  used  in  their 
code  (with  piecewise  bilinear  basis  functions)  converges  to  the  norm  of 
the  true  error.  We  have  not  yet  proven  convergence  for  our  error 
estimators. 

In  Chapter  3,  we  develop  several  heuristic  models  of  error  behavior 
which  lead  to  adaptive  refinement  procedures  (cf.  Brandt  [16],  Babulka 
and  Rheinboldt  [7]).  We  present  the  asymptotic  expected  error  and  work 
behavior  for  three  problem  classes  consisting  of  problems  with: 

1.  smooth  solutions. 

2.  solutions  with  point  singularities. 

3.  solutions  with  line  singularities. 

Except  in  the  presence  of  line  singularities,  the  expeoted  error 
converges  to  zero  faster  than  inverse  polynomially  with  respect  to  work. 
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In  Chapter  6,  we  discuss  some  of  the  computational  aspects  of  the 
methods  presented  in  Chapters  2-5.  Te  consider  efficient  assembly  of 
the  finite  element  system  (cf.  Bisenstat  and  Schnltz  [18],  Veiser, 
Eisenstat,  and  Schnltz  [31]),  and  efficient  construction  of  the  error 
estimators.  Ve  present  the  asymptotic  complexity  of  nested-disseetion- 
type  Gaussian  elimination  for  the  three  classes  of  problems  considered 
in  Chapter  5.  The  operation  counts  and  storage  for  the  problems  with 
singularities  are  linear  (optimal-order)  in  the  number  of  elements  N,  in 
contrast  to  the  operation  counts  and  storage  for  uniform  meshes,  which 
are  proportional  to  N3'2  and  N  log(N)  respectively  (e.g.,  George  [23]). 

In  Chapter  7,  we  present  numerical  results  obtained  using  prototype 
codes.  Ve  present  the  resulting  error  and  error  estimator  behavior  for 
problems  from  the  three  olusses  considered  in  Chapter  5.  The  error 
estimators  are  usually  accurate  to  within  a  factor  of  two.  In  some 
cases,  the  error  estimators  appear  to  converge  to  the  norm  of  the  true 
error  as  the  mesh  size  decreases. 

I  am  deeply  grateful  to  my  advisor.  Professor  Randolph  E.  Bank,  for 
suggesting  the  topic  of  this  research  and  working  closely  with  me  on  its 
realization.  I  am  honored  to  know  him  and  to  have  worked  with  him.  I 
would  like  to  thank  my  other  advisors.  Professors  Martin  H.  Schultz  and 
Stanley  C.  Eisenstat,  for  helping  to  guide  my  research,  and  for  their 
unfailing  intelligence,  integrity,  and  support. 

Many  other  past  and  present  members  of  the  numerical  analysis  group 
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at  Tale  have  contributed  much  to  my  enj  oyaent  of  graduate  achool  and  my 
appreciation  and  understanding  of  numerical  analysis.  They  are  (in 
order  of  appearance  at  Tale)  Andy  Sherman,  John  V.  Lewis,  Rati  Chandra, 
Rob  Schreiber,  Jack  Perry,  Dave  Fyfe,  Trond  Steihang,  Craig  Donglaa, 
Howie  Elman,  Ren  Jackson,  Tony  Chan,  John  Kindle,  and  Benren  Zhn. 

I  am  gratefnl  to  many  other  membera-at-large  in  the  numerical 
analysis  community.  I  am  especially  gratefnl  to  Professor  Nary 
F.  Wheeler,  for  her  encouragement  and  guidance  when  I  was  an 
undergraduate,  and  Professor  Ivo  Babutka,  for  his  excellent  and 
influential  work. 

I  am  grateful  to  the  Office  of  Naval  Research  (ONE  Grant 
N00014-76-C-0277.  FSU-ONR  Grant  F1.N00014-80-C-0076) ,  the  Department  of 
Energy  (DOE  Grant  DE-AC02-77ET53053) ,  the  National  Science  Foundation 
(Graduate  Fellowship),  and  Tale  University  for  their  financial  support. 

Most  of  all,  I  am  grateful  to  my  family,  and  my  loving  wife  Mary 

Ann. 


Loeal  auk  nflantit 


2.1  Iatroduetion 

la  this  chapter,  n  consider  bisection-type  local  aesh  refiaeaeat 
with  square  eleaents.  A  staple  aesh  of  this  type  is  shown  in  Figure 
2-1. 


Figure  2-1:  A  bisection-type  locally-refined  aesh 
This  type  of  aesh  refiaeaeat  has  been  intensively  investigated  in 


recent  years,  aainly  by  Babutks  and  Rheinboldt  and  their  students 
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(e.g.,  [4]).  The  data  structures  and  computational  algorithms  required 
to  implement  the  finite  element  method  with  general  meshes  of  thia  type 
are  fairly  complicated  (see  Rheinboldt  and  Mesrtenyi  [25]  and 
Gannon  [21]).  Thus,  Bank  and  Sherman  [14]  and  Van  Rosendale  [30] 
suggest  using  restricted  classes  of  meshes  which  can  be  implemented  with 
much  simpler  data  structures. 

In  Section  2.2,  we  present  some  standard  terminology  for  dealing 
with  bisection-type  local  mesh  refinement,  and  define  the  class  of 
"1-irregular"  meshes,  a  variant  of  a  class  of  meshes  suggested  by  Bank 
and  Sherman  [14].  In  Section  2.3,  we  consider  a  variant  of  a  mesh 
refinement  rule  suggested  by  Bank  and  Sherman  which  can  convert  any 
bisect ion-type  locally-refined  mesh  into  a  1-irregular  mesh  with  no  more 
than  a  constant  times  as  many  elements,  such  that  the  resulting  finite 
element  matrix  (with  bilinear  basis  functions)  has  no  more  than  a 
constant  number  of  nonseroes  in  any  row.  Conversely,  in  Section  2.4,  we 
show  that  there  are  simple  bisection-type  locally-refined  meshes  such 
that  the  resulting  finite  element  matrices  are  essentially  dense.  Te 
also  consider  alternative  mesh  refinement  rules.  In  Section  2.5,  we 
present  data  structure  details  for  handling  1-irregular  meshes. 

2.2  Terminology 

Ve  now  consider  bisection-type  local  mesh  refinement  with  square 
elements.  For  simplicity,  we  take  the  domain  to  be  the  unit  aquare, 

0  [0.11X10.1]  -  ((x,y):  Oixil,  Oiyil], 


with  boundary  dO.  The  square  region  E  ■  [xg.Xg+hgJXIyg.yg+hg]  it  called 
an  elenent.  The  elements 


!,B' W21XIyB+V2*  W  *  1 W2 ' W  X [yE+ V2 ' VV  ' 
lV  VV2]XIyE*yE+V21  *  [  V  V2 '  VV  X  tyE'yE+V21 

are  called  the  sons  of  E>  and  E  is  called  the  father  of  its  sons.  A 
mesh  N  is  defined  to  be  an  nnordered  set  of  elements.  An  element  in  N 
is  called  refined  if  its  four  sons  are  in  N.  and  unrefined  otherwise. 

The  class  of  bisection-type  locally-refined  meshes  with  square 
elements,  or  admissible  meshes  (Babulka  and  Bheinboldt  [8,  10]).  is 
recursively  defined  by  the  following  two  rules: 

1.  {0}  is  an  admissible  mesh. 

2.  If  a  mesh  N  is  an  admissible  mesh,  and  E  is  an  unrefined 
element  in  N,  then  the  mesh  {M.  the  sons  of  E}  is  an 
admissible  mesh. 

The  level  L  of  E  is  defined  to  be  the  number  of  generations  between 
Q  and  E.  Thus,  in  an  admissible  mesh  on  0,  h^  *  (1/2)L.  A  neiehbor  of 
E  is  defined  to  be  another  element  at  the  same  level  as  E  that  shares  a 
side  with  E. 

A  corner  of  an  unrefined  element  is  called  a  vertex.  A  vertex 
which  it  a  oorner  of  each  unrefined  element  it  touches  is  called  a 
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IinlU  vertex.  All  other  vertices  ere  celled  irr staler*.  The 
irregularity  index  (Babnlka  and  Rheinboldt  [10])  of  a  aesh  is  defined  to 
be  the  aaxiaua  nuaber  of  irregular  vertices  on  a  side  of  aa  oarefiaed 
eleaent.  Ve  call  a  aeah  with  irregularity  index  £  k  a  k— irreaular  aeah. 
A  (regular)  vertex  oa  90  is  oalled  a  boundary  vertex.  All  other 
vertices  are  called  interior. 

In  figures  depicting  aeshes, 

0  denotes  a  regular  vertex, 
o  denotes  a  vertex  (regular  •**’  irregular), 

•  denotes  an  irregular  vertex, 

(v)  at  a  vertex  denotes  function  value. 

For  exaaple.  Figure  2-1  depicts  an  adaissible  2-irregular  aesh  with  4 
refined  eleaents,  13  unrefined  eleaents  (one  labelled  E) ,  12  regular 
boundary  vertices,  3  regular  interior  vertices,  and  6  irregular  vertices 
(one  labelled  V) . 


This  is  the  naae  used  by  Bsbuhka  and  Rheinboldt  [10].  Bank  and 
Sheraan  [14]  use  the  naae  "green  vertex".  Gannon  [21]  uses  the  naae 
"inactive  vertex". 
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Let  M  be  an  admissible  mesh.  A  function  g  is  called  piecewise 
bilinear  on  M  if,  on  any  unrefined  element  E  in  N, 

g(x,y)  -  gCxg.yg)  (Xg+hg-xXyg+hg-y)  /  h| 

+  gtxg+bg.yg)  (x-xgjtyg+hg-y)  /  h£ 

+  «(VyE+V  (xE+Vx)(r"yE)  1  4 

+  8(xE+hE'yE+V  (*“*E,<y“yB)  1  4  * 


Let  the  finite  element  snace  S  (-  S(M))  be  the  space  of  continuous 
piecewise  bilinear  functions  on  M.  Let  the  regular  vertices  for  M  be 
V1,...,VN  .  Let  Bj(x.y)  (•=  B(Vj)(x,y))  be  the  function  in  S  which 
satisfies  the  Lagrange  conditions 


BI(V 


'IJ 


<  1  if  I  -  J 
0  if  I  *  J, 


for  J*1,...,N.  Bj  is  well-defined  using  induction  on  element  level. 
Moreover,  {B^}  is  a  basis  for  S,  since  each  Bj  is  in  S,  the  Bj's  are 
linearly  independent,  and  each  function  in  S  is  a  linear  combination  of 
the  Bj’i, 

The  sunnort  of  Bj,  supp(B^),  is  defined  to  be  the  union  of  all 
unrefined  elements  E  such  that  Bj  is  not  identically  xero  in  B.  In 
general,  the  support  of  Bj  may  be  non-oonvex  or  non-simply-connectsd. 
Alternately,  the  standard  tensor-product  basis  functions  (Bj)  may  be 
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defined  for  S,  so  that  each  supp(£(V))  is  a  rectangle.  Most  of  the 
considerations  applying  to  {Bj}  also  apply  to  {£j} .  However,  snpp(B(V>) 
is  always  contained  in  supp(£(V))  (see  Fignre  2-2),  and  we  know  of  no 
computational  advantage  for  (g^) .  Thus,  we  restrict  our  attention  to 


(0)-(0) - (0) 

I  I  I 

(0)-(0) (1/2)  I 

III  I 

(0)  (1/2)  (1) - (0) 

I  I  I 

I  I  I 

I  I  I 

(0) - (0) - (0) 

B(l/2.1/2) 


(0)-(0)-(0) - (0) 

III  I 

(0)(l/4)(l/2)  I 
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(0)  (1/2)  (1) - (0) 

I  I  I 

I  I  I 

I  I  I 

(0) - (0) - (0) 

£(1/2, 1/2) 


Figure  2-2:  Lagrange  versus  tensor-product  basis  functions 

2.3  The  1-irregular  rule  and  soae  of  its  eonsequenees 

Given  an  admissible  aiesh,  we  force  possible  additional  refinement 
by  applying,  as  many  times  as  possible,  the  following  vsrisnt  of  a  rule 
of  Bank  and  Sherman  [14] : 


refine  any  unrefined  element  with  more  than  one  (2.1) 

irregular  vertex  on  one  of  its  sides. 

Ve  call  rule  (2.1)  the  1-lrresnlar  rule.  The  1-irregular  rule  must 
terminate,  since  it  never  increases  the  maximum  element  level  in  the 
mesh.  After  the  1-irregular  rule  terminates,  the  resulting  mesh  is 
1-irregular.  Conversely,  if  any  refinement  forced  by  the  1-irregular 
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rale  is  omitted,  the  resulting  mesh  is  not  1-irregnler.  Thus,  the 
resnlting  mesh  hss  the  fewest  elements  of  sny  1— irregular  mesh 
containing  the  original  mesh. 

The  resulting  configurations  (up  to  symmetry)  of  the  support  of 
B(V)  for  interior  V  are  depicted  in  Figure  2-3.  Because  of  the 
1-irregular  rule,  further  refinement  of  an  element  in  one  of  cases  I-VI 
produces  another  of  cases  I-VI,  possibly  at  the  next  finer  level. 

Bach  unrefined  element  E  in  the  support  of  B(V)  satisfies  either 

A  :  V  is  a  corner  vertex  of  E, 
or 

B  :  V  is  not  a  corner  vertex  of  E,  but 

V  is  a  corner  vertex  of  the  father  of  E,  and 

V  is  next  to  an  irregular  corner  vertex  of  E. 

Thus,  it  is  simple  to  determine  whioh  basis  functions  are  nonsero 
in  E.  Figure  2-4  depicts  a  sample  E  in  its  father,  f(E),  and  numbers 
the  6  relevant  vertices.  B^  and  B^  are  nonsero  in  E  (V^  is  regular 
since  f(B)  is  refined).  Bj  (B^)  is  nonsero  in  E  if  V2  is  regular 
(irregular).  Bg  (B^)  is  nonsero  in  E  if  Vj  is  regular  (irregular). 

Iliere  are  exactly  four  basis  functions  nonsero  in  B.  Since  the 
element  stiffness  matrices  are  small  (4X4) ,  and  the  mappings  from  local 
to  global  basis  functions  are  simple,  the  resulting  finite  element 
system  is  straightforward  to  assemble. 
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Figure  2-3:  Coaf iguratioas  of  elements  in  snpp(B(V)) 
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Figure  2—4:  f(E) 

Furthermore,  the  basis  functions  are  locally  independent: 

in  each  unrefined  element  E,  (2.2) 

(Bjlg  :  E  in  supplB^)}  are  linearly  independent. 

where  *l£  denotes  restriction  to  element  E.  Ve  can  show  that  (2.2) 
holds  if  and  only  if 

no  irregular  vertex  touches  unrefined  elements  (2.3) 

at  three  different  levels. 

For  example,  in  the  mesh  shown  in  Figure  2-1.  irregular  vertex  V  touches 
unrefined  elements  at  levels  1.  2.  and  3,  and  the  restrictions  to 
element  E  of  the  five  basis  functions  nonzero  in  E  are  linearly 
dependent . 

Let  N(V)  be  the  number  of  basis  functions  with  support  intersecting 
supp(B(V)).  where  V  is  a  regular  vertex  in  a  1-irregular  mesh.  N(V)  is 
the  number  of  nonzeroes  in  the  row  of  the  finite  element  matrix 
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corresponding  to  B(V) .  The  minimum  possible  N(V)  for  as  interior  vertex 
is  5,  and  the  maximum  possible  N(V)  is  13  (Fignre  2-5).  N(V)  is  9  for 
interior  vertices  in  a  locally  uniform  mesh.  N(V)  for  boundary  vertices 
can  be  as  small  as  4.  N(V)  is  6  for  boundary  side  vertices  in  a  locally 
uniform  mesh. 
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Figure  2-5:  Minimum  and  maximum  N(V)  for  interior  vertices 

Many  refinement  sequences  arising  in  practice  naturally  yield 
1-irregular  meshes,  so  the  1-irregular  rule  introduces  so  additional 
refinement.  Ve  now  show  that,  given  an  arbitrary  admissible  mesh,  and 
applying  the  1-irregular  rule  as  many  times  as  possible, 

the  number  of  elements  forced  to  refine  by  the  (2.4) 

1-irregular  rule  is  no  more  than  eight  times  the  number 
of  refined  elements  in  the  original  mesh. 

Bence  the  number  of  basis  functions  introduced  by  the  1-irregular  rule 
is  no  more  than  a  constant  times  the  number  of  original  basis  functions. 

Suppose  there  are  LMAX  levels  in  the  mesh.  Let 
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0(L)  :=  {original  refined  elements  at  level  L} , 

R(L)  :=  {elements  at  level  L  refined  by  the  1-irregnlar  rule}. 

It  suffices  to  shov  that  each  element  E  in  R(L)  shares  a  corner  vertex 
with  at  least  one  element  F  in  O(L) .  We  do  this  using  induction  from 
Lr*LMAX  back  down  to  0.  There  are  clearly  no  elements  in  R(LMAX) . 

Suppose  the  statement  holds  for  all  levels  with  index  greater  than  L. 
Since  E  is  in  R(L),  the  1-irregular  rule  forces  E  to  refine  becanse  G,  a 
neighbor  of  E,  has  a  son,  H,  which  is  refined  at  level  L+l,  and  shares  a 
side  with  E.  If  H  is  in  0(L+1),  then  G»F  is  in  O(L) ,  and  (2.4)  holds. 
Otherwise,  H  is  in  R(L+1).  By  induction,  there  is  an  element  I  in 
0(L+1)  which  shares  a  corner  vertex  with  H.  Then  the  father  F  of  I  is 
in  O(L)  and  shares  a  corner  vertex  with  E  (Figure  2-6).  So.  (2.4)  holds 
in  all  cases. 
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Figure  2-6:  H  in  R(L+1) 


We  do  not  know  if  (2.4)  is  sharp.  The  largest  ratio  of  enforced 
refinement  to  original  refinement  we  have  encountered  is  5.75. 
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In  summary,  given  any  admissible  mesh,  the  1-irregular  rule  forces 
additional  refinement  such  that 

the  resulting  1-irregular  mesh  has  no  more  than  a  (2.S) 

constant  times  the  original  number  of  basis  functions, 

and 

the  resulting  element  stiffness  matrices  are  of  (2.6) 

bounded  size  (4X4),  and  simple  to  assemble, 

and 

the  number  of  nonzeroes  in  any  row  of  the  resulting  (2.7) 

finite  element  system  is  bounded  independent  of  the  mesh. 

2.4  Alternative  mesh  refinement  rules 

In  this  section,  we  show  that  a  refinement  rule  like  the 
1-irregular  rule  is  needed  for  properties  (2.6)  and  (2.7).  We  also 
consider  alternative  mesh  refinement  rules. 

Suppose,  because  of  a  sharp  point  singularity,  that  a  mesh  results 
from  refining  only  elements  which  contain  the  point  P  -  (2/3, 2/3). 

Since 

2/3=1-  1/2  +  1/4  -  1/8  +  1/16  -  .  .  .  , 
the  mesh  is  formed  by  successively  refining  the  alternately  upper  right 
or  lower  left  son  of  the  smallest  refined  element.  If  there  are 
LMAX  >>  1  levels  in  the  mesh,  the  active  vertices  are  the  ten  boundary 
vertices  and  the  LMAX  interior  vertices 

\  -  <  PL  •  Pl  >  •  PL  *  1  1-0  M . 
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Let  element  E  be  the  smallest  element  containing  P.  For  any  ,  it  can 
be  shown  by  indnetion  on  I  that  B(V^)  is  nonzero  at  the  points  (Pj_j,Pj) 
and  if  and  only  if  I  >  L.  Setting  I  -  LMAI,  since 

*PLMAX-1'PLMAX^  an^  ^PLMAX'PLMAX-1^  ar*  oorner*  °*  E*  E  mnst  be  in  the 
support  of  each  B(V^) .  Thus,  the  element  stiffness  matrix  for  E,  and 

the  resulting  finite  element  system,  are  essentially  dense.  Since  the 
mesh  is  2-irregular,  this  example  disproves  the  conjecture  (Babulka  and 
Miller  [4],  page  17)  that  (2.7)  holds  for  k-irregular  meshes  with  k 
bounded. 

Assembly  of  the  finite  element  sytem  is  fairly  complicated  (e.g., 
Gannon  [21],  Rheinboldt  and  Mesztenyi  [25]).  If  Gaussian  elimination  is 
used  to  solve  the  system,  the  solution  time  is  asymptotically 
proportional  to  LMAI3 .  By  contrast,  if  the  1- irregular  rule  is  applied, 
assembly  is  simple.  If  the  vertices  are  ordered  according  to  the 
maximum  levels  of  the  elements  in  their  supports,  solution  time  is 
asymptotically  proportional  to  LMAI  (see  Section  6.3). 

In  the  previous  example,  (2.6)  and  (2.7)  can  be  insured  without 
enforcing  additional  refinement,  by  using  modified  basis  functions  which 
are  identically  zero  in  the  element  containing  the  singularity.  Figure 
2-8  depicts  the  modified  basis  funotion  g ( ) .  When  there  are  two  sharp 
point  singularities  at 

P1  -  (9/16  +  (2/ 3 ) / 16 ,  10/16  +  (2/3 ) / 16 ) 
and 

P2  -  (9/16  +  (2/3)/16 ,  13/16  +  (2/3)/16), 


however 
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Figure  2-7:  Singularity  at  P  =  (2/3, 2/3),  LMAX  *  4 
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Figure  2-8:  A  modified  basis  function 
,  the  support  of  the  basis  function  for  must  asymptotically 


intersect  the  supports  of  roughly  half  the  basis  functions  in  S,  since 
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it  can  only  be  zero  at  one  point  along  the  line  segaent  froa  (1/2, 1/2) 
to  (1/2,1).  Thna,  with  any  rnle  for  foraing  baaia  functions,  it  aay  be 
advantageous  (purely  froa  aatrix  sparsity  considerations)  to  force 
additional  refinement. 

Ve  briefly  aention  alternatives  to  the  1-irregular  rnle.  Another 
rule  which  insures  properties  (2.3),  (2.6),  and  (2.7)  (when  applied  as 
many  tines  as  possible)  is: 

refine  all  unrefined  side-  and  corner-  neighbors  of  aay  (2.8) 

eleaent  which  has  grandsons  (sons  of  a  son). 

Rnle  (2.8)  is  equivalent  to  a  condition  in  Van  Rosendale  [30],  page  59. 
With  this  rule,  only  cutes  I-V  in  Figure  2-3  are  possible,  and  (2.5), 
(2.6),  and  (2.7)  hold  for  both  [Bj)  and  {Jlj} .  However,  rule  (2.8)  aay 
force  aore  additional  refineaent,  and  henoe  require  aore  work,  than  the 
1- irregular  rule,  and  is  no  easier  to  iapleaent. 

Conversely,  any  rule  which  aay  force  less  additional  refineaent 
than  the  1-irregular  rule  does  so  only  when  neighboring  eleaents  with 
size  ratios  at  least  4:1  occur.  For  ezaaple,  if  we  apply  the  rule 

if  an  irregular  vertex  touohes  unrefined  eleaents  at  three  (2.9) 
different  levels,  refine  the  eleaent  at  the  aiddle  level, 

as  aany  tines  as  possible  (refining  eleaents  with  lowest  level  first), 
the  resulting  aeah  will  satisfy  (2.6),  sad  uz  have  fewer  eleaents  than 
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the  resulting  1-irregulsr  mesh.  Unfortunately,  (2.S)  end  (2.7)  any  not 
hold:  for  a  problen  with  a  sharp  point  singularity  at  (1/2+e ,1/2+e) ,  as 
e  approaches  zero,  the  nnaber  of  nonzeroes  in  the  row  of  the  finite 
eleaent  system  for  B(l/2,l/2)  may  be  arbitrarily  large. 

Figure  2-9  depicts  the  results  of  applying  the  1-irregular  rule, 
rule  (2.8),  and  rule  (2.9)  to  the  mesh  in  Figure  2-7. 

2.5  A  aesh  data  structure 

The  data  structure  we  use  to  represent  1-irregular  aeshes  is  a 
variant  of  a  data  structure  presented  in  Bank  and  Sheraan  [14].  The 
data  structure  has  two  aain  coaponents:  a  refined-element  tree  IRCT, 
stored  in  a  4XNIRCT  array,  and  a  vertex  list  IVERT,  stored  in  a  2XNIVERT 
array.  For  each  refined  eleaent  E,  there  is  a  4X4  node  in  IRCT 
containing  pointers  to  the  father  of  E,  to  any  of  the  four  sons  of  E 
which  are  refined,  and  to  the  nine  vertices  which  are  corners  of  the 
sons  of  E.  For  each  interior  vertex  V  which  bisects  an  eleaent  side, 
there  is  a  2X1  node  in  IVERT  containing  pointers  to  the  elements  with 
sides  bisected  by  V.  Other  types  of  vertices  have  other  inforaation 
stored  in  IVERT.  Figure  2-10  depicts  a  sample  refined  eleaent  E,  its 
node  in  IRCT,  and  a  sample  node  in  IVERT. 

The  neighbors  of  any  eleaent  can  be  easily  found.  For  instance,  in 
Figure  2-10,  the  pointer  location  for  the  bottoa  neighbor  of  Sj  is  in 
the  node  for  E,  while  the  pointer  location  for  the  right  neighbor  of  Sj 
is  in  the  node  for  the  right  neighbor  E'  of  E.  The  node  for  B’  (if  it 
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Figure  2-10:  Nodes  in  IRCT  end  IVERT 
exists)  is  pointed  to  in  the  node  for  Wj  in  IVERT. 

Moreover,  snppose  the  entire  mesh  is  rotated  90,  180,  or  270 
degrees.  Then  IVERT  and  the  top  row  of  IRCT  remain  the  same,  and  the 
last  three  rows  of  each  node  in  IRCT  are  rotated.  For  example,  if  the 
mesh  is  rotated  90  degrees  counterclockwise,  the  node  for  E  becomes 


f  (E) 


S 


1 


Thus,  the  neighbors  of  any  son  of  E  can  be  found  with  the  same  procedure 
used  to  find  the  neighbors  of  Sj  ,  as  long  as  indices  into  the  last 
three  rows  of  IRCT  are  incremented  by  the  appropriate  value  mod  4. 
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This  aesh  data  (true tare  sake*  it  relatively  eaay  to  construct  the 
finite  eleaent  system  for  S,  to  evaluate  the  appxoxiaate  solution,  and 
to  refine  the  aesh.  More  data  structure  considerations  are  in  Bank  and 
Sheraan  [14] . 


1  re  - - 


CHAPTER  3 


Local  order  refinement 


3.1  Introduction 

Traditional  finite  element  methods  for  elliptic  problems  esqploy 
finite  element  spaces  with  loir-order  (e.g.,  bilinear)  basis  fnnetions. 
These  spaces  are  appropriate  when  the  solutions  to  the  elliptic  problems 
are  not  smooth.  When  the  solutions  are  smooth,  standard  approximation 
theory  indicates  that  high-order  piecewise  polynomial  basis  fnnetions 
can  yield  much  faster  convergence.  Thus,  it  is  advantageous  to  be  able 
to  use  both  low-order  and  high-order  basis  functions. 

Local-order  finite  element  spaces  have  received  recent  attention  by 
Babulka,  Katz,  and  Szabo  [3,  11,  12]  and  Babulka  and  Dorr  [2].  Full 
local-order,  loeal-mesh  methods  in  more  than  one  dimension  have  not  yet 
been  considered. 

Ia  this  chapter,  we  extend  the  finite  element  spaces  with 
1-irregular  meshes  considered  in  Chapter  2  to  include  basis  functions 
with  locally-varying  polynomial  orders.  In  Section  3.2,  we  oonsider  the 
case,  and  in  Section  3.3,  we  briefly  consider  spaces  of  smoother 
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functions.  In  Section  3.4,  we  present  staple  data  structures  for  the 
spaces  considered. 

3.2  The  Lagrange  ease 

Let  M  be  a  given  1-irregular  nesh.  In  this  section,  we  construct  a 
space  S  (=  S(M))  of  C®  piecewise  polynosiials  such  that,  in  each 
unrefined  eleaent  E  in  M, 

saooth  functions  are  locally  approxinated  to  order  kg  2  2.  (3.1) 

Our  local  order  refineaent  follows  the  basic  approach  of 
Taylor  [29]  and  Babukka,  Katz,  and  Sxabo  [3].  On  [0,1],  let  the 
polynoaials  (p^x) :  i“l». . .  ,kaax]  satisfy 

P^x)  =  1-x, 

P2(x)  =  x, 

and  for  i  >  2, 

p^(x)  is  a  polynoaial  of  order  exactly  i, 

P4(0)  *  P1<1)  “  0. 

Sinoe  [p^,...,p^]  is  a  basis  for  the  polynoaials  of  order  k  >  1,  the  p^ 
are  called  (C°)  hierarohlc  polynoaials. 

E  E 

On  eleaent  E,  p^(x)  and  Pa<y)  ere  defined  by  napping  [0,1]  onto 
lWV  aad  tyE'yE+hE1 : 
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Djp^(x)  :=  Djpi((x-x£)/hE)  hE  j. 

and 

D^pm(y)  :=  D^PB((y“yE)^hE)  **£ 
j=0 , . . . , kmax-1 . 

Let  p,.(t)  be  the  Legendre  polynomial  on  [0,1]  of  order  i.  Since 
Pj(t)  =  1,  and 

fl  pi  pj  =  0  •  1  *  i‘ 

a  convenient  choice  for  {p ^ }  is 

pi(*>  =  /q  P^tt)  dt  ,  i  2  2.  (3.2) 

We  define  S  to  be  the  span  of  three  types  of  basis  functions 
(Bj(x.y)}:  vertex,  side .  and  element . 

The  vertex  basis  functions  are  the  piecewise  bilinear  basis 
functions  from  Chapter  2,  with  Bj(Vj)  **  8jj  for  all  regular  vertices  Vj. 

The  side  basis  functions  are  products  of  piecewise  linear  functions 
in  one  direction  with  quadratic  or  higher-order  hierarchic  polynomials 
in  the  other  direction.  Each  side  s  of  an  unrefined  element,  not 
strictly  contained  in  a  side  of  another  unrefined  element,  has  side 
basis  functions  with  hierarchic  polynomials  of  order  up  to 
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(3.3) 

kg  : =  max {k^  :  s  contains  a  side  of  unrefined  eleaent  E) . 

£ 

A  basis  function  for  side  s  in  Figure  3-1  is  h(x)pj(y). 


Figure  3-1 x  A  side  basis  function 

The  eleaent  basis  functions  for  an  unrefined  eleaent  E  are  the 
functions 

p^(x)  p^(y)  :  i»a  >  2,  i+a-1  ±,  kg  .  (3.4) 

Eleaent  basis  functions  for  E  are  nonzero  only  in  E,  and  only  occur  when 
kg25. 

Since  the  aesh  is  1-irregular,  each  basis  function  Bj  for  S  is  of 
Ex  Ev 

the  fora  p,  (x)  p  y(y)  in  each  unrefined  eleaent  E,  where  Ex  is  either  E 

i  a 

or  f(E),  and  Ey  is  either  E  or  f(E).  Moreover,  (3.3)  and  (3.4)  iaply 
that  in  each  eleaent  E,  any  aonoaial  x*ya  with  i+a-1  i  kg  is  a  linear 
combination  of  restrictions  to  E  of  basis  functions,  so  (3.1)  holds. 


*  •  '  ■'  -X 


Figure  3-2  depicts  the  values  of  the  basis  functions  in  an  element 
in  s  sample  local-order  finite  element.  Note  that  many  of  the  basis 
functions  are  of  order  greater  than  kg  =  3. 

3.3  The  Hermite  case 

The  construction  in  the  last  section  extends  immediately  to  spaces 
of  globally  C°  1  functions  for  any  integer  a  2.  The  hierarchic 
polynomials  satisfy  Hermite  interpolation  conditions  on  their  first  a 
derivatives  at  0  and  1.  so  kg  2.  2a.  Vertex,  side,  and  element  basis 
functions  are  defined  as  in  the  case.  There  are 

2 

a  vertex  basis  functions  for  each  regular  vertex, 

(kg-4a) (kE-4a+l)/2  element  basis  functions  for  element  E 
(none  if  kg  <.  4a) , 

and 

(k#-2a+l-i)  side  basis  functions  for  side  s, 
and  condition  (3.1)  holds  as  before. 

3.4  A.  mesh  data  strueture  for  looal  order  refinement 

In  this  section,  we  extend  the  mesh  data  structure  from  Section  2.5 
to  handle  local-order  finite  element  spaces. 

For  each  unrefined  element  E,  we  store  -kg  in  IRCT,  in  the  location 
in  the  node  for  f(E)  whioh  will  oontain  a  pointer  to  the  node  for  E  if  E 


is  refined 
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mesh  : 
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basis  functions  in  element  E  : 


Figure  3-2:  Looal-order  basis  functions 


Let  the  function  MBV  nap  tha  baais  functions  for  S  to  the  regular 
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vertices  in  the  nesh,  such  that 

1.  MBV  maps  each  vertex  basis  function  to  its  vertex. 

2.  MBV  naps  each  side  basis  function  for  side  s  to  the 
last-created  regular  vertex  at  an  end  of  s. 

3.  MBV  maps  each  element  basis  function  for  element  E  to  the 
center  vertex  of  f(E) . 

A  sample  mesh  indicating  MBV  is  shown  in  Figure  3-3. 

V  :  vertex  basis  functions 
S  :  side  basis  functions 
E  :  element  basis  functions 

@ :  basis  functions  mapped 
to  vertex  V 


Figure  3-3:  Basis  functions  mapped  to  vertices 

Ve  order  the  vertices  and  basis  functions  for  S  so  that  the  basis 
functions  mapped  by  MBV  to  vertex  Vjy  are 

BMVB(IV)  '  BMVB(IV)+1  '  •"  *  BMVB(IV+1)-1  ’ 
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and  we  associate  a  number  IBFT(IBF)  with  each  basis  function  , 

encoding  the  spatial  relation  to  the  vertex  MBV ( Bjgp)  (e.g.,  left,  lower 
right,  etc.). 

With  the  data  structure  from  Section  2.S,  and  two  additional  arrays 
containing  MVB  and  IBFT,  we  can  easily  determine  the  basis  functions 
which  are  nonzero  in  a  given  unrefined  element. 

As  in  Chapter  2,  this  mesh  data  structure  makes  it  relatively  easy 
to  construct  the  finite  element  system  for  S,  to  evalnate  the 
approximate  solution,  and  to  refine  the  mesh. 


CHAPTER  4 


A  posteriori  error  estimators 


4.1  Introduction 

The  traditional  error  estimates  for  the  finite  element  method  for 
elliptic  problems  are  a  priori  (e.g.,  Ciarlet  [17]).  predicting  the 
expected  rate  of  convergence  of  the  error  to  zero,  but  not  saying  much 
about  the  actual  error  for  a  given  finite  element  space. 

In  recent  years,  loeally-computable  a  posteriori  error  estimators 
have  been  developed,  mainly  by  Babutka  and  Rheinboldt,  which  estimate 
the  actual  error  for  a  given  space.  The  landmark  paper  in  this  area  is 
by  Babutka  and  Rheinboldt  [8],  giving  a  method  of  constructing  a 
locally— computable  error  estimator  which  is  within  a  multiplicative 
constant  of  the  norm  of  the  actual  error.  Further  results  for 
one-dimensional  problems  are  given  by  Babuhka  and  Rheinboldt  [6,  7,  5] 
and  Reinhardt  [24].  Further  results  for  two-dimensional  problems  are 
given  by  Babulka  and  Miller  [4]:  under  suitable  assumptions,  they  prove 
that  the  error  estimator  used  in  their  code  (with  piecewise  bilinear 
basis  functions)  converges  to  the  norm  of  the  actual  error  as  the  mesh 
size  goes  to  zero. 
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In  Section  4.2,  we  present  some  terminology  for  discussing 
s  posteriori  error  estimators  for  two-dimensional  Neumann  problems.  In 
Section  4.3,  we  briefly  describe  the  error  estimator  presented  in  [8]. 


In  Section  4.4,  we  present  a  new  error  estimator  which  has  two  main 
advantages  over  the  estimator  in  [8] : 

1.  Under  reasonable  assumptions,  the  new  error  estimator  is 
asymptotically  an  upper  bound  on  the  norm  of  the  true  error. 

2.  The  new  error  estimator  can  be  oomputed  locally  in  each 
element.  (The  one  in  [8]  must  be  computed  over  more 
complicated  local  regions.) 

In  Section  4.3,  we  prove  the  upper-bound  property  of  the  new  error 
estimator,  and,  under  suitable  assumptions,  we  show  that  the  new  error 
estimator  is  no  larger  than  a  multiplicative  constant  times  the  norm  of 
the  actual  error.  In  Section  4.6,  we  present  several  cheaper 
alternative  error  estimators.  In  Section  4.7,  we  extend  the  estimators 
to  handle  more  general  boundary  conditions  and  differential  operators. 
In  section  4.8,  we  summarize  the  results  in  this  chapter. 

In  Chapter  7,  we  present  numerical  evidence  that,  in  some  cases, 
our  error  estimators  converge  to  the  norm  of  the  actual  error  as  the 
mesh  size  goes  to  zero.  We  have  not  been  able  to  prove  this 
convergence.  Currently,  we  can  only  prove  convergence  for  our  error 
estimators  for  one-dimensional  problems  [151. 


4  ■> '«» •  *-  %• 
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4.2  Terminology 

Consider  the  model  Neumann  problem 


Ln  :=  -  Dx(a(x,y)Dxn)  -  Dy(a(x,y)Dyu)  +  b(x,y)u  -  f(x,y)  (4.1) 


in  the  interior  of  Q, 
du/dn  =0  on  dQ, 

where  ■  and  d/dn  denote  derivatives  in  x,  y,  and  the  outward 

normal  to  0  respectively,  a(x.y)  and  b(x.y)  are  in  L  (0) ,  f(x.y)  is  in 
L2(Q),  0  <  a(x,y)  i  a.  and  0  <  Jt  i  b(x,y)  £  b.  Variants  of  (4.1) 
are  discussed  in  Section  4.7. 


Let  u  be  an  open  subset  of  0  with  piecewise  smooth  boundary  du. 

00 

For  a  non-negative  integer  k,  and  a  function  v  in  C  (w) ,  we  define  the 
norm 


MvMk.W  :=  1  i“0  }  J-0  fl  Jv(X'y»2  d*  dy* 

Let  H  («)  denote  the  completion  of  C  («)  with  respect  to  the  norm 

11*11.  .  H*(w)  is  equivalent  to 

x,« 

(v  in  L2(«)  :  llvll.  <  •} , 


if  derivatives  are  defined  in  the  distributional  sense  (e.g.. 

Ciarlet  [17],  page  114).  Let  Hq(»)  denote  the  completion  of  the 
00 

functions  in  C  («)  with  compact  support  in  w  with  respect  to  the  norm 

11*11.  When  k  is  omitted,  k  ■  0  is  assumed.  For  v,w  in  H^(«) 

x,« 
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(=  L^(a>) ) ,  let 

(v,w)w  :*  /£  v(x,y)  w(x,y)  dx  dy. 

For  v,w  in  H*(«),  let 

e  (v,w)  :*=  (aD  v,D  w)  +  (aD  v,D  w)  +  (bv,w) 

(t)  x  X  tt  y  y  m  (0 

If  w  is  the  interior  of  Q,  a  (u,y)  is  obtained  by  Multiplying  (4.1) 

u 

function  v,  and  integrating  over  m  using  Green's  theorea.  Let  the 
2 

energy  norm  1 1 1 v 1 1 1  :*  a  (v,v) .  By  the  assumptions  on  a(x,y)  and 

O 

b(x,y) ,  there  exists  C  >  0  independent  of  v.  such  that 

c"1  IMIj,.  i  III.III,  it  IItII1># 

for  all  v  in  H1 («) . 

2 

On  the  boundary  3u>,  for  v  and  w  in  L  (d*i),  lot 

<v’w>a«  L  y(o)  wU)  do* 

where  do  is  arc  length.  We  define  the  norm 


C  denotes  a  generic  constant,  not  necessarily  the  same  in  eaoh 


by  a 


(4.2) 


instance. 
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,v,a«  :=  <v'v>a«  • 


Let  d/dn  denote  differentiation  in  the  direction  of  the  outward  noraal 
to  (ii,  with  the  limit  in  the  differentiation  taken  over  points  inside  w. 


For  brevity,  we  sometimes  replace  a  by  its  closure  in  naaes  for 

spaces,  norms,  and  inner  products:  for  ezaaple,  1 1 v 1 1  „  stands  for 

a,  c 

1 1 v 1 1 ^  interior(E)  .  When  the  subscript  u  is  omitted  in  noras  and  inner 
products,  a  =  interior(O)  is  assuaed. 

Given  a  1-irregular  aesh  N  with  unrefined  eleaents  {£} ,  we  choose  a 
finite  element  subspace  S  (=  S(M>)  of  H*(Q)  as  in  Chapters  2  and  3.  For 
simplicity,  we  assuae  S  is  based  on  C®  hierarchic  polynoaials,  as  in 
Section  3.2.  We  also  choose  a  larger  finite  element  space  S  (=  S(M)) 
containing  S,  constructed  in  the  sane  manner  as  S,  such  that  in  each 
unrefined  element  E, 

Since  S  contains  S,  and  the  basis  functions  {B^}  are  hierarchic,  the 
basis 


smooth  functions  are  locally  approximated  to  order  kg 


{Bj  :  Bj  in  S] 
for  S  contains  the  basis 


{Bj  :  Bj  in  S} 
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for  S.  For  each  unrefined  element  E,  we  define  S_  to  be  the  space  with 

£> 

basis  functions 

{Bjlg  :  E  in  suppfB^,  Bj  in  S} , 

where  tBjIg  :  E  in  snpp(B^)}  are  linearly  independent  because  of  (2.2), 
(3.3),  and  (3.4),  and  we  define  Sg  to  be  the  space  with  basis  functions 

(Bj.|g  :  E  in  supp(Bf),  Bj  in  SJ . 

Ve  define  Sg-Sg  to  be  the  span  of  the  functions 

(Bjlg  :  E  in  suppCB^),  Bj  in  S,  Bj  not  in  S) . 

For  any  function 

V  “  5  BjaS  YI  BI*E 
in  Sg  ,  we  define 

IE<v)  ^  B^sS  Wj  Bjlg  (4.3) 

in  Sg  .  Then  wlg(w)  is  in  S^-Sg  ,  and  for  any  function 

V  "  I  BjeS  yI  BI 
in  S,  the  function 

I<T)  ='  l  Bj.S  TI  BI 
is  in  S,  and  is  equal  to  Ig(v)  in  E. 


t 
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Let  N(E)  denote  the  set  of  unrefined  elements  which  intersect  with 
E  at  more  than  one  point.  For  example,  in  Figure  3-2,  N(E)  is  the  set  of 
labelled  elements. 

For  convenience  in  discussing  functions  in  the  product  space 
:=  {v:  v  in  H1 (E)  for  all  E} , 
let 

a(v,w)  :=  ^  £  ag(v,w)  for  v  and  w  in  HM(Q) , 

and 

lllvlll2  :«  ^  E  lllvlllg  for  v  in  B^(Q) . 

The  weak  form  of  (4.1), 

a(u.v)  =  (f,v)  for  all  ▼  in  H*(Q), 

has  a  unique  solution  u  in  H*(0)  (e.g.,  Ciarlet  [17],  page  19).  Ve 
impose  the  weak  form  of  (4.1)  in  S  to  get  a  finite  element  approximation 
D  in  S,  where 

a(U,v)  =  (f,v)  for  all  ▼  in  S. 

The  error  for  S  is  e  :*  u-D.  Similarly,  U  in  S  is  required  to  satisfy 
a(D,v)  ■  (f,v)  for  all  v  in  S, 
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and  the  corresponding  error  for  S  is  e  :*  n-U. 

In  this  chapter,  we  consider  ways  of  estimating  the  non  of  the 
error  lllelll.  Let 

e  :*  0-0, 

Ve  as some  the  satnration  condition  that  there  exists  0  £  y  <  1  snch  that 

lllelll  1  y  lllelll.  (4.4) 

If  n  is  sufficiently  smooth,  then  y  goes  to  zero  as  h  goes  to  zero. 

Thus,  since 

lllslll  1  lllelll  i  (l-?)-1  HUM. 

we  concentrate  our  attention  on  estimating  lllslll. 

4.3  A  Diriohlet  a  poateriori  error  eatimator 

Traditionally,  the  error  estimates  for  the  finite  element  method 
have  been  a  priori,  predicting  the  expected  rate  of  convergence  of 
lllelll  to  zero  (nsnally  as  h  :*  max(hg)  goes  to  zero),  but  not  saying 
much  about  the  actual  error  for  a  given  S. 

In  [8] ,  Babulka  and  Rheinboldt  present  a  posteriori  error 
estimators  constructed  basieally  as  follows.  Let  (Bj)  be  a  subset  of 
the  basis  functions  for  S  which  form  a  partition  of  unity  on  0.  Por 
each  J,  let  Qj  supp(IO,  and  let  tj  be  the  solution  in  of  the 


local  Diriohlet  problem 


Coder  suitable  assumptions,  Babulka  and  Rbeinboldt  show  that  there  exist 


Cj  >  0  and  >  0  independent  of  S  such  that 

C1  nielli2  <  }  J  IIUjIll^  i.  C2  lllelll2. 

In  practice,  since  H^(Qj)  is  infinite-dimensional,  e^  is  not 
computable.  The  most  straightforward  computable  approximation  of  e^.  is 
cj  in 

Sj  :=  (p(x,y)  :  p  in  S,  supp(p)  contained  in  Qj) 
defined  by  the  equations 

Bfl  ^eJ*v^  *  <f’v>0  "  -q  for  a11  v  in  Sj  .  (4.5) 


Then 

C1  HlelN2i}j  i  C2  lllelll2  .  (4.6) 

^  2  1/2 

The  quantity  (  )  j  1 1 1 e j 1 1 1 Q  )  is  called  an  error  estimator. 

J 

Each  «j  in  (4.5)  is  a  projection  of  e  in  the  energy  norm  onto  a  set 
of  functions,  so  the  right  hand  inequality  in  (4.6)  arises  naturally, 
while  the  left  hand  inequality  is  harder  to  prove. 


4.4  A  N«omb  a  posteriori  error  estimator 


The  motivation  for  a  Neumann  a  posteriori  error  estimator  is  that 
it  is  more  Vmportant  to  have  an  npper  bound  on  lllslll  than  a  lover 
bound.  Ve  construct  an  error  simnlator  s,  such  that  e  is  the  projection 
of  s  onto  S  with  respect  to  the  energy  norm.  Then  Nielli  j,  lllslll 
arises  naturally.  Following  Babulka  and  Rheinboldt  [10],  ve  call 
lllslll  an  error  estimator,  and  for  each  element  E,  ve  call  1 1  lei llg  an 
error  indicator. 

To  motivate  our  construction  of  e,  suppose  u  is  in  H2 (E) .  By 
Green's  theorem  on  element  E, 


Sg(e,v)  »  F£'v)  :■  (f.v>E  +  <*du/dn,v> -  Sg(D,v)  (4.7) 

for  all  v  in  H* (E) .  For  a  given  v,  ve  can  compute  every  quantity  in 
(4.7)  except  adu/dn  on  the  boundary  of  E.  Suppose  E  shares  side  s  vith 
element  E' .  On  s,  a  possible  approximation  to  adu/dn|g  is  the  average 

where  ng  is  the  outvard  normal  from  E. 

Sinoe 


value  of  adU/dng  on  E  and  E' , 


nE,s  "  "  V1.  ' 


ve  define  the  approximation 


[adU/dnJ  lB  (ill  {  sdD/dnlg  -  adU/dnlR,  }  ,  s  not  in  dll  (4.8) 


'E' 
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Note 


[»3D/a»]lE  -  -  [*ao/aa]lE,  . 

The  resulting  approximation  to  Fg(v)  i* 

Fg(v)  :*  (f,v)g  +  <ta3D/dn] ,v>gE  -  aE(C,v).  (4.9) 

Since  ag(*,*)  is  elliptic,  we  can  define  an  error  simulator  for  e  in  E 
by  replacing  Fg  with  Fg  on  the  right  hand  side  of  (4.7).  However,  as  hg 
approaches  zero,  the  finite  element  matrices  for  the  error  simulator 
approach  an  indefinite  matrix  corresponding  to  the  case  when  b(x,y)  =  0. 

Thus,  we  formulate  error  simulator  equations  which  are 
automatically  consistent  when  b(x,y)  -  0  and  aE(*,*)  is  indefinite.  We 
approximate  a  by  the  function  e£  in  Sg  defined  by  the  equations 

ag(eE.v)  =  Fg(v-IE(v))  for  all  v  in  Sg  ,  (4.10) 

where  Ig(v)  i*  defined  by  (4.3).  Since  8g(»,*)  is  elliptic,  8g  exists 
uniquely.  Since  ▼  need  not  vanish  on  3E,  (4.10)  is  a  Neumann  problem. 
Moreover,  since  llE  is  in  Sg  , 

IE(1)  -  1. 

so  the  consistency  condition  £g(*g>l)  *  0  is  automatically  satisfied. 
Equations  (4.10)  are  equivalent  to  choosing  s_  in  S_  suoh  that 
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aE  ®E'V>  "  |  FE(v)  for  v  ln  SE'SE 
JO  for  v  in  Sg 

The  corresponding  matrix  equation  is  A  v  *  f  for  the  coefficient  vector 
V  of  «E  ,  where 


00  a 

A-  1 

*  I 

*2  1 


0  I 
I 


the  coefficients  of  the  basis  functions  for  S  are  in  v,  ,  and  the 

E  1 

coefficients  of  the  basis  functions  for  Sg-S^  are  in  , 

We  define  the  value  of  the  error  simulator  s  at  (x,y)  in  0  to  be 
^(x.y)  S’gtx.y), 

where  E  is  an  unrefined  element  containing  (x,y) .  When  there  is  no 
danger  of  confusion,  we  refer  to  «g  as  e.  Note  that  «  is  in 
not  necessarily  in  1^(0). 

4.5  Properties  of  the  Neumann  error  estimator 
The  most  important  property  of  c  is  that 

i»l:  III. Ill  1  1117111. 

Proof :  Suppose  v  is  in  S.  Since  [a9U/9n]|g  -  -  [a9U/9n]|g,  ,  and 
[adU/dn]  *  0  on  9Q, 


2  g  Fg(v>  ■  2  e  ^f«v^E  ~  Sg(U,v)}  =  (f#v)  -  a(U, v) . 
Since  I(v)  (page  38)  is  in  S, 

5  E  V*E<v))  "  1  E  V1(v))  "  "  «(®»X(t))  -  0. 


Thus, 


a(s,v)  *  a(e,v)  =  (f,v)  -  a(D,v)  “  ^  E  ^(v-IgCv))  “  ■(«»▼)» 
and,  in  particular, 

1 1 1 e 1 1 ! 2  =  a(e.s)  =  a(s,e)  £  lllelll  lllelll. 


By  (4.4),  we  have 

Corollary  4.2;  (1-y)  lllelll  £  lllelll. 

This  is  a  alobal  inequality.  Ve  cannot  generally  expect 
llle|||g£  C  lll.lllg  :  suppose  u  is  saooth,  S  contains  bilinears, 

hE“ 

a  finite  eleaent  approxiaation  to  the  solution  to 


h  in  region  R  in  0,  and  h^  *  h  >>  h  away  froa  R.  In  region  R, 


.11) 


.12) 


U  is 


L(u^)  “  f(x,y)  in  the  interior  of  R 


(4.13) 
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Let  E  be  an  element  in  the  interior  of  R.  Since  the  error  simulator 
equations  in  E  are  local,  they  are  the  aaae  for  (4.13)  aa  for  (4.1). 
Since  (Theorem  4.6)  lllelllg  is  bounded  by  local  quantitiea  behaving 
like  IllU-Uglllg  .  and  the  boundary  error  in  (4.13)  polluter  the 
solution  throughout  R.  ve  expect 

lllclllE  »  lllelllg  . 

(In  practice,  this  may  be  acceptable,  sinoe  the  error  ia  bigger  away 
from  R  due  to  the  coarse  mesh,  and  the  small  size  of  the  error  simulator 
is  in  accord  with  the  fact  that  no  more  mesh  refinement  should  occur.) 

By  the  triangle  inequality  and  (4.12), 

lllslll  -  III. Ill  1  1 1 1 s-a 1 1 1  -  inf (I  I le-f  I 1 1 :  f  in  S) . 

so  we  can  get  an  a  posteriori  upper  bound  for  1 1  Is  1 1 l-l 1 1.1 1 1  by 
measuring  the  difference  in  the  energy  norm  between  a  and  any  function 
in  S.  For  instance,  by  averaging  a,  we  can  construct  an  approximation 

sea 

a  to  a  in  S:  s-s  approximates  the  sice  of  the  jumps  in  a  across 
ave  ave 

element  boundaries.  If 

1 1 1  a-a  II  lr  Nielli 
ave 

with  y  <  1,  then 

lllslll  i  (1-y)"1  lllslll. 

Of  course,  the  jumps  in  s  can  be  made  as  small  as  desired  by  adding 


penalty  tern*  to  equations  (4.10),  bat  the  resulting  equations  are 
non-local . 

We  now  prove  a  local  a  priori  upper  bound  for  llleglllg  in  terms  of 
quantities  depending  on  e.  We  assume  there  exists  a  constant  kmax  such 
that 

kg  <_  kmax  independent  of  S.  (4.14) 

We  start  with  a  few  lessees. 

LfiSfi&lti:  |I|IE(JE)I|IE  i  C  lllsglllg  . 

Proof:  Suppose  v  is  in  Sg  .  Since  the  dimension  of  Sg  is  finite  and 
is  linear, 

(bIE(v).IE(v))E  S.  Cj  (bv,v)g  , 

where  is  independent  of  v  and  hg  .  (Otherwise,  there  would  exist  w 
in  Sg  such  that  (bw,w)g  =  0  and  (blg(w) ,Ig(w) )g  #  0,  which  leads  to  a 

is  finite,  Ig  is 

linear,  and  Ig(l)  =  1. 

(aDllE(v).DxIE(v»E  +  <‘Ve<T>'Ve(v))E  * 

C2  (  (aDiv,Dxv)g  +  (aDyv,Dyv)g  ), 
where  C.  is  independent  of  v  and  h_  .  Thus, 


contradiction.)  Similarly,  since  the  dimension  of  Sg 
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llllE(v)|||g  i  ninfCj.Cj)  llUlllg  . 

Let  v  =  8g  . 

Lemma  4.4:  I Sg-Igtfig) lflE  i  C  hg^2| 1 1 Bgl 1 1£  . 

Proof:  The  trace  inequality 

*8E~IE(8E)  -  C  thE  ^eE-IE(eE)  **E  +  ^  *  *  8E_IE(8E>  ^1,E}  (4*15) 

follows  directly  fron  Theorem  3.10  in  Agaon  [1]. 

Suppose  v  is  in  Sg  .  For  0  i  x,y  i  1,  let 

v0(x.y)  :=  v(xE  +  hE  x.  yg+hjf), 

yx.y}  :=  IE(v)(xE  +  \  x,  yE  +  ^  y>  - 
Since  the  dimension  of  SE  is  finite.  Ig  is  linear,  and  Ig(l)  "  1, 

*,to'i(iiiq  i  c  *l*DiT«**i  *  *llV,a**a  ’■ 

and  thus 

I  lvwIE<y)  I  |g  IChJ  IMI^g  • 

where  C  is  independent  of  ▼  and  hg  .  The  result  follows  by  letting 
▼  =  s_  in  (4.15)  and  using  (4.2)  and  Leans  4.3. 

fi 

Leaaa  4.5:  If  n  is  in  H^E), 


-  W  -  • 
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I-/**!,,  i  c  S,'1'2  <  "•ll2.K  >,/2' 


Proof :  Since  a  is  in  H^(E),  e  is  in  H^(E).  The  result  follows  directly 
from  Theorem  3.10  in  Agmon  [1]. 


The  following  theorem  shows  that  the  error  is  large  in  local 
regions  where  the  error  indicators  are  large  (and  n  is  sufficiently 
smooth) . 


Theorem  4.6:  If  n  is  in  (E')  for  all  E'  in  N(E) ,  then 
llleglllg  1  C  ^  b»bN(E)  *  +  hE*  ^*^2,B'  1 


Proof:  For  any  ▼  in  Sg  , 

*E*“e'v*  “  ^f,v“IE^v^E  +  »▼“!£<▼)  >gB  ~  ag(D,w-Ig(v)) 

*  Fg('r-Ig(v))  -  <adn/0n*y“lg(y)>ag  +  <[adU/3n] ,w-Ig(y)  >aE 
■  ag(e,y-Ig(v))  -  <[a9e/3n] ,y-Ig(v)>aE  . 

Let  y  ■  «g  .  By  the  Cauchy-Sohwarx  inequality  and  the  definition  of 

M. 

II  l«E,,,E2  i  l,l»H,E  ,,,*e'IE(‘e),ME  + 

1  5  E'eN(E)  *ta</anl3E'  }  * ‘b’W  *dE  * 

By  Lenaa  4.3  and  Leaaa  4.4, 
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Proof:  It  follows  fro*  Theorem  3.10  is  Agmon  [1]  tad  (4.14)  that 

i  c  hjf1  lll.llj  g  ♦  4  ll.ll*#I  ) 

Then,  using  sde/dn  “  ade/dn  +  ads/dn  in  Equation  (4.16), 

1 1 1  a  ^  1 1 1  g  C  lllalllg  +  C  hg  5  E’sN(E)  + 

C  5  E’sN(E)  *I,**1,B'  * 

Squaring,  s naming  oyer  E,  using  the  triangle  inequality,  and  noting  that 
the  maximum  number  of  N(*)*s  in  vhioh  any  element  appears  is  bounded, 

lllelll  i  C  lllelll  ♦  C  I  lei  Ij  +  C  (JE  hg  ladi/dnljg  )1/2. 

The  Theorem  follows  using  (4.2),  noting  that  lllelll  £  lllelll,  and 
using  (4.17). 

4.6  Alternative  error  eatlmatore 

In  praotice,  there  are  several  alternative  error  estimators  that 
are  less  oostly  to  ooavute  than  lllelll. 

4 >6*1  Ignoring  irregular  earners 

Ve  oaa  compute  an  error  simulator  ia  element  B  by  replaeimg  Ig(v) 
with  an  analogous  function  I*(v)  which  interpolates  v  loeally  in  E.  The 
resulting  error  simulator,  ia  effect,  ignores  the  irregularity  of  the 
corners  of  B,  and  its  system  is  slightly  cheaper  to  assemble  than  the 
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system  A  v  *  f.  In  practice,  the  resulting  error  estimator  appears  to 
be  at  least  as  accurate  as  I  list  II.  However,  Equation  (4.11)  in  the 
proof  of  Theorem  4.1  break  down  when  v  is  a  side  basis  function  for  a 
side  s  containing  an  irregular  vertex  (Figure  4-1):  I  (v),  the  analogue 
to  I(v),  is  not  continuous  across  s.  Thus,  we  do  not  consider  this 
alternative  error  estimator  further. 


(0)- 

1 

- (0) - 

1 

- (0) 

1 

(0)- 

1 

- (0) - 

1 

- (0) 

1 

1 

(1)1 

1 

1 

(1)1 

1 

(0)- 

- (1) 

1 

(0)- 

- (0) 

1 

1 

(1)1 

1 

1 

(1)1 

1 

1 

1 

1 

1 

1 

1 

(0)- 

- (0) - 

- (0) 

(0)- 

— (0) - 

- (0) 

v  I*(v) 

Figure  4-1:  Ignoring  the  irregular  oorners 

4.6.2  Using  the  matrix  for  a  local  Poisson  problem 

We  can  form  and  solve  error  simulator  equations  using  a  matrix 
corresponding  to  an  associated  local  Poisson  problem.  For  v  and  w  in 
H*(E),  let 

Ag(v,w)  :*  g(E)  {  (D^v.D^g  +  ^DyV,Dyw*E 

IJ.lvlJ.lg  5“  lg(v.v), 

where  *(E)  inf(a(x,y)  in  E} . 

Instead  of  choosing  8g  in  Sg  by  (4.10),  we  can  choose  ig  in  Sg  such 


that 


S3 


Vse-’1  '  "f5 

L 


E(y)  for  all  v  in  Sg-S 


for  all  ▼  in  Sg  , 


and 


(4.18) 


(8g.l)g  -  0.  (4.19) 

(Ve  inpose  (4.19)  because  £g(l,l)  “  0.)  We  construct  the  error 
sianlator  e  and  tbe  error  indicators  (I  Mai  Mg)  from  in  the  saaie  way 
we  constructed  e  and  (lllslllg)  from  eg  . 

Equations  (4.18)  correspond  to  the  natrix  equation  A  v  *  f  for  the 
coefficient  vector  v  of  8g(x,y),  where 


A 


K  I 
’  I 

*2  1 


v 


I 

I 

I 


f 


I  0  I 

I  I  . 

I  f,  I 


the  coefficients  of  the  basis  functions  for  Sg  are  in  ,  and  the 
coefficients  of  the  basis  functions  for  Sg-Sg  are  in  Vj  .  Since 
IjJiljUg  "  lila+ClHg  for  any  constant  C,  we  oan  set 

*1,J  “  ^J*l  “  8J.l  * 


reordering  the  unknowns  if  neoessary.  so  that  A  i*  positive  definite. 

For  fixed  kg  ,  A  depends  only  on  the  regularity  of  the  corners  of  B,  and 
not  on  hg  . 
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Theorem  4.8 :  If  a(x,y)  is  Lips chi tz-cont inuous ,  then 

lllel  HE  1  lll«lllE  <  (l+Odig))  Mlelllg 

as  hg  goes  to  zero. 

Proof:  By  (4.10)  and  (4.18).  since  ±g(v,v)  £  Bg(v.v)  for  all  ▼  in  S. 

IlKlIlJ  -  J^(a.e)  i  MleMlg  INllNg  i  MleMlg  lllllllg  . 

Similarly. 

MlaNIg  =  aE(a.s)  =  a^J.a)  i  lllelllg  lllalllg  . 

and  if  a(x,y)  is  Lipschi tz-cont innona  in  E,  so  that 
a(x.y)  i  *(E)  +  O(hg)  in  E.  then  by  (4.19). 

M-ME  i  c  hE  MiMl,E  ' 


lllalllg  <  (l+0(hg))  MlalMg 
as  ‘.'g  goes  to  zero. 

Since  mesh  refinement  occurs  in  a  very  regular  fashion,  we  can 
precompute  the  factorization  of  and  only  incur  the  cost  of  a 
backsolve  plus  computation  of  the  right  hand  side  for  each  element, 
the  hierarchic  polynomials  (p^J  are  chosen  by  (3.2),  then  since 


If 


0  ,  i  4  j  ,  max( i, j)  >  2 


ss 


A  is  sparse.  In  practice,  we  replaoe  the  quantity  *(E)  with  the  wore 
convenient  quantity  a(Xg+hg/2,yE+hg/2) .  Then 

lllslllg  i  (1+0(1^))  MIsMIg  i  (l-KXhg))  lll?lllE. 


4.5.3  Using  a  smaller  matrix 

We  can  also  compute  8g  in  Sg-Sg  by  the  equations 

%<®E'v)  “  ^e(v)  for  a11  v  in  ®E_SE  '  (4.20) 

constructing  the  error  simulator  a  and  the  error  indicators  (I  I  lei  1 1 _} 

—  -  £ 

from  ep  in  the  same  way  we  conatruoted  e  and  (lllelll  }  from  e„  . 

Equations  (4.20)  correspond  to  the  matrix  equation  A  v  «  f  ,  or 

£  £  £ 

A  v  =  f,  for  the  coefficient  vector  v  of  i  ,  where 

£ 


I  A  0  I  I  0  | 

A  I  1  I  .  v  :«  I  | 

I  «  *2  I  I  \  1 


f 


I  o  I 

I  I  . 

I  f,  I 


We  now  show 

Theorem  4.9:  There  exists  a  constant  C  independent  of  hg  such  that 

lil«MIEI  C  MUINg  . 


Proof :  The  left  hand  inequality  follows  from 

MUIMg  -  vtav  -  vTf  -  vtAj  -  igU.e)  i  IMsMIg  MlaMlg  . 
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The  following  approach  to  the  right  hand  inequality  is  dne  to 

Eisenstat  [19].  Since  A  is  positive  definite,  nsing  block  Ganssian 

elimination,  A.  v„  =  f„  ,  where 
4  2.  2 

*4  <A2  -  *3  V 

is  symmetric  and  positive  definite.  Then 

iiuiuj  -  f 2  a;t  f2  . 

Similarly,  A^  is  positive  definite,  and 

Ul«l He  =  4  A2T  f2  * 

Since  A.  and  A.  are  f ini te-dimens ional , 

4  2 

-  2  T  -T  T  -T 

C  *  :=  max  v1  A/  v  /  v1  A  1  v  < 

v*0  4  2 

In  particular,  let  v  =  f2  . 

For  example,  if  E  has  no  irregular  corners,  a(E)  =  1,  S,„  contains 
bilinears,  and  Sg-Sg  contains  the  biquadratic  basis  functions 

( Wx)  <yE+Vy)  '  4  * 

(*“*E>  (yE+hE"y)  (ywyE)  '  4  ' 

(WX)  (X"XE)  (yE+Vy>  '  4  * 

( WX)  (x-xE)  (y-yE)  '  4  ' 
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then 


1 

6 

0 

0 

0 

1 

1  0 

0 

0 

0 

1 

1 

0 

4 

-1 

-2 

l 

1  1 

1 

-1 

-1 

1 

0 

-1 

4 

-1 

1  /  6  . 

A,  -  I  -1 

1 

1 

-1 

1  /  12  » 

1  1 

0 

-2 

-1 

4 

I 

3  1  -1 

-1 

1 

1 

1 

1  13 

0 

2 

0 

1 

1  0 

13 

0 

2 

I 

A,  “  1  2 

0 

13 

0 

1  /  90  , 

^  1  0 

2 

0 

13 

1 

:  =  (ii/6) 

1/2  ~ 

1.354. 

4.6.4  Using  the  rea ideals  and  jnps  in  normal  derivatives 
We  may  choose  not  to  solve  any  linear  systems. 

By  Green's  theorem,  for  v  in  Sg-Sg  , 

Sg(U,v)  *  (LU, v)g  +  <adU/9n,v>3E  , 

so 

~E^8E’v^  *  ^ f-LU,v)g  +  < [a3U/9n]-e3U/dn,v>9g  . 

In  particular, 

IMsgl  Mg  -  (f-LU,iE)E  +  <[a«U/an]-adU/an,;E>aB  . 


and  so 
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MIbeI!Ie  <  CjOig.^tE))  1 1  f-LDl  Ig  + 

CjO^.atE))  I  [a3U/dn]-adU/dnlgE  , 


(4.21) 


where 


^(hg.aCE))  :=  sup  I  lv|  Ig/M  IvIHg  » 


V8S+ 


Cjdig.AtE))  :=  sup  IvIgg/IMvMlg  , 
V8S+ 


and  S+  :=  {v  :  v  =  vQ(x,y)+C,  vQ  in  Sg-Sg  ,  v  *  0,  and  (v,l)E  **  0} , 


For  example,  if  E  has  no  irregular  corners,  Sg  contains  bilinears, 
and  Sg-Sg  contains  the  biquadratic  basis  functions  on  page  56,  then  S+ 
is  the  span  of 

**E  *XE+hE_X*  *yE+llE-y*  "  1^12» 

hg  3  (x-Xg)  (yE+hg-y)  (y-yg)  -  1/12, 

hg  3  (xg+hg-x)  (x-xg)  (yg+hg-y)  -  1/12, 

^  3  (xE+bE-x*  (x-xE^  (y"yE*  "  1^12* 


and 


C^hg.itE))  -  hg  /  (22  A(E))1/2, 
C2(hg,A(E))  -  {  3  hg  /  (11  A(E)>  }1/2. 


Babulka  and  Rheinboldt  [9]  suggest  error  indicators  when  S  contains 
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piecewise  bilinear  baaia  functions  which  mount  to  approximating 

lllelllg  -  h£/(24A(E))  ||f-LU||2+  (4.22) 

hE/(6i(E))  l[aaU/an]-adD/3nlgE 

for  interior  elements.  Since  1/22  -  1.1  X  1/24  and  3/11  «*  1.6  X  1/6. 
the  two  schemes  are  in  rough  agreement. 

If  the  mesh  is  uniform.  S  contains  piecewise  bilinear  basia 
functions.  S  contains  piecewise  biquadratic  basis  functions,  u  is 
sufficiently  smooth,  and  h  goes  to  zero,  then  combining  Corollary  4.2, 
Theorem  4.8,  Theorem  4.9,  and  (4.21), 

1 1  lei  1 12  i  (I-})"2  J  E  ((124(E))"1/2hEllf-Lll|lE  + 

(2jl<E))-1/2hE1^2|  [ adU/dn] -a8U/ dn I  2 , 

where  y  goes  to  zero.  Similar  results  hold  for  higher-order  elements 
and  nonuniform  grids. 

In  numerical  tests,  we  have  obaerved  that  quantities  of  the  form 
J  B  C1(hE.A(E))  Hf-LUll|  +  C2(hB,i(E))  I  [adU/dn] -adU/dn I 2e 
or 

J  E  {^(hg.^E))  Ilf-LUllg  +  C2(hB,i(E))  I  [  adD/dn]  -adU/dn  I  #E) 2 

2 

only  seem  to  approach  I  llel  II  to  within  a  multiplicative  constant 
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depending  on  u.  In  Chapter  7,  we  present  numerical  results  using  the 
error  estimator 


i-  ‘5 


2  ,1/2 
e„  )  , 


E  =E 


where 


:=  |(kE){h|/s(E)  Ilf-LUllg  +  4hE/A(E)  I ta0O/an]-adU/0nl ,(4.23) 


snd 


ftkg)  :=  (k^kg+lHkg+l))"1.  (4.24) 

e„  is  chosen  to  extend  formula  (4.22).  and  to  perform  well  on  the 
problem 

2  2 

-Du-Du*f(x,y)  in  0  ,  u  =  0  on  80, 

x  y 

and  hE  are  uniform,  and  h  goes  to  zero. 

4.7  Extensions 

If  b(x,y)  =  0  in  Q,  in  order  to  insure  nonsingular  linear  systems, 
we  impose  the  Neumann  consistency  conditions  (u,l)  **  0,  (0,1)  *  0,  and 
( 8 ,1 ) E  ■  0.  We  still  have 

lllvlllg  i  C  llvllj  E  for  all  v  in  H*(E), 

and  globally, 

Mel  I ^  1  C  III. III. 


when  u  =  sin(nx) sin(ny) ,  kg 
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Since  <IE»1)E  “  0  M*EME  i  C  hg  ll?Ell1<E  , 

ii«ii1  Bic  iii;iiie 

for  sufficiently  small.  (Similarly,  llellj  E  i  C  1 1  lei  I lg  .)  Thna, 
Theorem  4.6  and  Theorem  4.7  hold  for  h  anff iciently  email. 

If  aE(*,*)  ia  of  the  form 

a  (v,w)  -  (a(x,y)  D  v,D  »)„  +  (a(x,y)  D  y,D  w)_  + 
n  x  x  e$  y  y  b 

(b(x,y)  v,w)g  +  (c(x,y)  D^y.wjg  +  (d(x,y)  Dyv,w>E  , 

then,  under  anitable  assumptions,  llsll*  g  £  C  aE(e,e)  and 

II  a  1 1 2  E  £  c  »*)  as  h  goes  to  xero  (see  Sohatx  [27]).  Thna,  Theorem 

4.6  and  Theorem  4.7  hold  for  fa  sufficiently  email. 

For  problems  with  non-homogeneous  Neumann  boundary  conditions 

adu/dn  ”  g(x,y)  on  90  ,  g  in  H^(9Q), 

if  ve  set 

a(u,v)  -  ( f ,▼)  +  <g,T>  for  all  y  in  H*(Q), 
a(0,y)  “  (f,y)  +  <g,y>  for  all  y  in  S, 

and 

[adU/dn]  :■  g(x,y)  if  a  is  on  9Q, 


then  we  obtain  the  same  results  as  before. 

For  problems  with  homogeneous  Dirichlet  boundary  conditions 
u  =  0  on  80, 

the  solution  u  in  Hq(Q)  satisfies 

a(u,v)  =  (f,v)  for  all  v  in  hJ(0). 

•W*  — 

Similarly,  the  equations  for  U,  s,  e,  and  e  are  now  posed  on  subspaces 
of  H^(Q) .  Since  all  the  functions  of  interest  now  vanish  on  dQ,  the 
other  results  follow  as  before.  Problems  with  non-homogeneous  Dirichlet 
boundary  conditions  which  are  exactly  satisfied  are  handled  similarly. 

When  Dirichlet  boundary  conditions  are  not  exactly  satisfied  by  D, 

— 

in  element  E  touching  dQ,  we  define  e^  in  S£  by  the  equations 

*E(Vv)  "=  ^  *or  ***  v  ^E  w,lich  o® 

Bp  ”  e  on  dQ. 

We  define  e  from  e„  as  before.  Let  e.  be  any  function  in  S  such  that 

K  D 

Equation  (4.19)  is  no  longer  needed  in  elements  touching  dQ. 


<3 
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goes  to  zero.  Since  e  is  nonzero  on  9Q,  e  should  be  in  Sg  rather  the 
Sg-Sv  when  E  touches  dO. 

Analogous  results  hold  for  e  and  s. 

As  in  Babulka  and  Rheinboldt  [10],  on  an  eleaent  side  s  of  E  which 
touches  90,  we  replace  the  ter® 

|(kE)  dl^/aCE)  l[a9U/9n]-a9U/9n|2 

with  the  tera 

3a(E)/h„  lei2 
IS  8 

in  foraula  (4.23). 


4.8 


The  error  indicators  discussed  in  this  chapter  are: 


llUglUg  ,  where  8g  is  in  Sg  and 


a„(8„.v)  =  /f. 


E  E 


.g(v)  for  all  ▼  in  Sg-Sg 


for  all  v  in  Sg  , 


IJ.liglJ.lg  *  w>>ere  *E  is  in  SR  ,  (Sg,l)g  -  0,  and 


E  '  -E'E 


hE(iE'v) 


Fg(v) 


1° 


for  all  ▼  in  Sg-Sg 
for  all  ▼  in  Sg  , 
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MUglNg  ,  where  t£  ia  in  Sg-Sg  and 

AE*eE'V*  *  £or  #11  V  ia  ®E-SE  ' 

and  &g  ,  where 

eg  -  |(kg)  (hJ/A(E)  ||f-LD||g  +  4hg/A(E)  I [adU/dn]-adU/dnl£g  }. 


The  corresponding  error  e«tl»ators  are 

iiiein  -  ( J  E  m;Em l  )in. 

MUNI  -  (  }  E  MUglNg  )1/2. 

MUNI  -  (  }  E  MUgMIg  )in. 


and 


«2.* 


2  .1/2 
E  }  ' 


e 

Te  define  the  corresponding  relative  errors  in  the  estiaiators 


Babnlka  and  Rheinboldt  call  (the  relative  error  +  1)  the  "efficiency 


index"  [10]  or  the  "effectivity  index"  [9]. 
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P  :=  (  IMelll  -  lllelli  )  /  Nielli, 

£  (  MUM  -  lllelli  )  /  lllelli, 

l  (  MleMI  -  lllelli  )  /  lllelli. 


end 


g  :«  <  e  -  lllelli  )  /  lllelli. 


Ve  have  shown,  for  sufficiently  smooth  n,  that  as  h  goes  to  zero. 


(1-;)  lllelli  1  MUM  1  C  lllelli. 

<lU  lllelli  i  lllelli  i  C  lllelli, 

<!-?>  C_1  lllelli  i  lllelli  iC  lllelli. 


and 


C  MUM  1  e. 


where  the  positive  constants  C,  C,  and  C  are  independent  of  h,  and  if 
is  sufficiently  smooth,  y  goes  to  zero  if  h  goes  to  zero. 


*.  - 


mat 


Heuristic  refinement  criteria 


5.1  Introduction 

The  galding  principle  in  choosing  s  finite  element  space  is  to 
minimize  the  norm  of  the  error  as  a  function  of  work.  However,  there  is 
not  usually  enough  information  available  to  choose  an  appropriate  space 
a  priori.  Instead,  it  is  necessary  to  choose  a  sequence  of  spaces 
adantivelv.  Each  succeeding  space  in  the  sequence  is  chosen  on  the 
basis  of  criteria  computable  from  the  preceding  space,  with  the  criteria 
tending  to  satisfy  the  guiding  principle  when  a  heuristic  model  of  error 
behavior  holds. 

When  the  polynomial  orders  of  the  basis  functions  are  uniform,  the 
resulting  refinement  algorithms  tend  to  approximately  equalize  the  sizes 
of  the  errors  in  the  elements  in  the  next  space.  In  finite  element 
methods,  for  example,  this  heuristic  refinement  criterion  is  used  in 
Fried  and  Tang  [20],  Babulka  and  Rheinboldt  [8],  and  Bank  and 
Sherman  [13].  When  the  polynomial  orders  of  the  basis  functions  vary 
locally,  the  heuristic  refinement  criteria  must  be  somewhat  more 
complicated  (e.g.,  Brandt  [16]). 
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In  Sections  5.2  -  5.5,  we  develop  two  heuristic  models  of  error 
behavior  for  local— mesh,  local— order  finite  element  spaces,  and  an 
adaptive  refinement  algorithm  which  uses  the  models.  In  Section  5.6,  we 
discuss  drawbacks  inherent  in  the  models.  In  Section  5.7,  we  present 
the  asymptotic  expected  error  and  work  behavior  for  three  problem 
classes  consisting  of  problems  with: 

1.  smooth  solutions. 

2.  solutions  with  point  singularities. 

3.  solutions  with  line  singularities. 

Except  in  the  presence  of  line  singularities,  with  appropriate  local 
mesh  and  local  order  refinement,  the  expected  error  converges  to  zero 
faster  than  inverse  polynomially  with  respect  to  work. 

5.2  A  model  of  work 

The  guiding  principle  in  choosing  a  "good"  finite  element  space  S 
can  be  stated  as 

minimize  I  Mali  I  as  a  function  of  V,  (5.1) 

where  1 1  loll  I  is  the  error  in  the  energy  norm,  and  T  is  the  work. 
However,  there  is  not  usually  enough  information  available  to  choose  an 
appropriate  space  a  priori.  Thus,  we  must  choose  a  sequence  of  nested 
spaces 

S,  contained  in  S.  contained  in  ...  S  , 

\  l  n 

with  finite  element  solutions  U, ,...,D  ,  taking  S  ■  S  (for  instance, 

X  n  n 

because  the  work  limit  is  reached).  We  choose  each  adantivelv.  on 

the  basis  of  information  obtained  from  S ^  . 


Since  the  spaces  are 
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nested,  the  errors  { eA}  and  work  counts  {W^J  satisfy 
III 111  2  llle2ll|  1  ...  2  IlleJII 
and 

^  ^  ^  ...  ^  . 

The  work  counts  can  represent  either  storage  or  coaputer  tiae.  If 

the  work  counts  represents  storage,  ¥  -  W  .  If  the  work  counts 

n 

represents  tiae,  ¥  -  ^  and  if 

Wi  ^  y  \+l  *  0  i  T  <  1  .  1  i  i  i  n,  (5.2) 

then  W  i  (1-y)  *  W  .  Thus,  ¥  is  usually  proportional  to  W  .  Since  it 
o  & 

is  possible  that  n  =  i+1,  in  the  reaainder  of  this  chapter,  we  take  as 
the  "local  guiding  principle"  a  local  version  of  (5.1): 

given  SA  ,  choose  Si+1  to  miniaize  IlleJII  I  (5.3) 

as  a  function  of 

(subject  to  the  possible  enforceaent  of  (5.2)). 

Let  S  be  the  current  finite  eleaent  space  in  the  sequence,  with 
aesh  N,  unrefined  eleaents  (E) ,  and  error  e.  Let  z  :■  (x,y)  denote  a 
point  in  0.  Let  h(x),  p(z),  and  w(p(z))  be  piecewise  constant  functions 
defined  by 
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h(z)  fey  =  the  local  mesh  size, 

p(z)  :=  pE  := 

w(p(z))  :=  »E  =  the  local  work  count, 

if  z  is  in  element  E.  We  make  the  heuristic  assumption  that  the  work  W 
for  S  is 

W  =  5  E  *E  ‘  (5,4) 

This  is  reasonable  if  W  is  the  dimension  of  S,  or  (since  the  mesh  is 
1-irregular)  the  storage  for  the  stiffness  matrix  for  S,  or  the  number 
of  arithmetic  operations  needed  to  solve  the  finite  element  system  for  S 
when  the  work  is  proportional  to  the  number  of  elements  in  the  mesh  (see 
Section  6.3).  It  is  not  reasonable  if  the  work  grows  faster  than  the 
number  of  elements:  then,  W  depends  on  S  globally  as  well  as  locally. 

5.3  A  model  of  predicted  error  behavior 

In  this  section,  we  construct  a  model  of  predicted  error  behavior 
based  on  standard  a  priori  error  estimates,  and  derive  heuristic 
refinement  criteria  which  tend  to  satisfy  (5.3)  when  the  model  holds. 

Let  m(z)  and  p(z)  be  piecewise  constant  functions  defined  by 

m(z)  :=  max{  k  :  Hull.  „  <  •  )  , 

K,  C 


kg-1  =  the  local  polynomial  degree. 


p(z)  :»  min(p(z) ,m(z)-l) 
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if  z  is  in  element  E.  m  indicates  the  maximum  local  smoothness  of  u  in 

a# 

E,  and  p  indicates  the  expected  convergence  rate  in  the  energy  norm 
nsing  polynomials  of  degree  p  (order  k^) .  Standard  a  priori  error 
estimates  (e.g.,  Ciarlet  [17],  Theorem  3.1.5  and  proofs  of  Theorems 
3.2.1  and  3.2.2)  yield 

a* 

1 1 le| 1 12  i  C(p(z) )  a(z)  G( z.p(z) )  h(z)2p*z^  dz,  (5.5) 

where  a(z)  is  the  coefficient  in  (4.1),  and 


G(  z 


}  J 


P+1  (  Dj  DP+1-j  u(z)  )2 


*  y 


is  in  L  (Q)  by  the  definition  of  p. 


Following  Brandt  [16],  Chapter  8,  we  make  the  crncial  assumption 
that,  for  heuristic  purposes,  (5.5)  holds  in  each  element,  with 
bflSAli.ty: 

M 

lllelll2  =  :*  C(p(z))  a(z)  G(z,p(z))  h(z)2p  **  dz,  (5.6) 

so  that 

Hum2  = «  =-JE«E . 

Ve  also  assume,  for  heuristic  purposes,  that 


a  and  G  are  approximately  constant  in  each  E 


(5.7) 
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C,  G,  and  W  are  (formally)  differentiable  in  h  and  p.  (5.8) 

Assumptions  (5.4),  (5.6),  (5.7),  and  (5.8)  comprise  onr  first  heuristic 
model  of  predicted  error  behavior.  The  form  of  (5.6)  is  important,  but 
not  the  specific  values  of  C(p)  and  G(z,p).  Ve  make  assumption  (5.8) 
only  for  convenience:  all  derivatives  with  respect  to  h  and  p  can  be 
replaced  by  divided  differences. 

In  the  following  paragraphs,  using  (5.4),  (5.6),  (5.7),  and  (5.8), 
we  estimate  the  marginal  efficiency  of  refinement  in  each  element  E, 
obtaining 

kjh)  =  -  (35/ahg)  /  OW/ahg)  (5.9) 

for  mesh  refinement  (h-ref inement) ,  and 

AgP)  Z  -  0{/3pE)  /  OW/apE)  (5.10) 

for  order  refinement  (p-ref inement) .  Let 
:=  max{AEh* ,\EP>) . 
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We  use  these  quantities  in  the  following  adaptive  refinenent  algorithm: 

1.  Find  the  unrefined  element  B  with  the  largest  Xg  . 

2.  Determine  the  refinement  cutoff  X^^  <  Xg  (e.g.,  as  in 
Section  5 .5) . 

3.  For  each  unrefined  element  E  with  Xg  2  X^^.  : 

a.  If  4b)  2  4P)  •  r«tln»  the  mesh  in  E. 

b.  If  Xgh)  <  XgP\  refine  the  order  in  E. 

4.  Take  S^+1  to  be  the  resulting  finite  element  space. 

We  now  determine  X^h>  and  X^p) .  First,  we  estimate  £  in  each 
unrefined  element  E.  Suppose  the  order  in  E  has  never  been  refined,  and 
(from  a  previous  finite  element  spaoe  in  the  sequence)  we  have  an 
approximation  to  the  square  of  the  error  in  f(E)  before  f(E)  was 

refined.  By  (5.6)  and  (5.7), 

*  2  «  ,2(p+l) 

0 

-  ln((g)  )  /  ln(4)  -  1. 

On  the  other  hand,  suppose  the  order  in  E  has  been  refined,  and  we  have 
an  approximation  to  the  square  of  the  error  in  element  B  with  kg  one 
less  than  its  current  value.  If 


’f(E)  *E 


* 


so  we  estimate 


,  -  (  U(«f(B)) 
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^  *E  ^  PCUT  (5.11) 

for  some  constant  p^  <  1  (e.g.,  as  in  Section  5.4),  there  may  still  be 
higher  derivatives  of  a  in  E  which  can  be  exploited,  so  we  estimate 
P  =  Pg  .  Otherwise,  we  retain  our  previous  estimate  of  p.  In  order  for 
our  estimate  of  p  to  be  logically  consistent,  we  restrict  it  to  the 
interval  (0,p£] • 

Using  our  estimate  of  p,  by  (5.7)  and  (5.9),  we  estimate 
4h)  =  -  (3(h|  a  C  G  h2p)/3h}  /  0(h£  w  h_2)/3h}  (5.12) 

=  -  {2  p  (a  C  G  h2p)  h_1)  /  {(-2)  w  h"3) 


=  (  P  /  »  ) 


in  element  E. 


There  are  two  cases  to  consider  in  estimating  X^P^ .  If  p  2  m-1, 
refining  the  order  will  not  change  the  error  at  all  (in  (5.6)).  Thus, 
4P>  =  0.  On  the  other  hand,  if  p  <  m-1,  by  (5.7)  and  (5.10), 

4P)  -  -  (a(«E)/ap)  / 

where  w'  : •  dw/dp.  If  u(z)  *  C  sin(Ox)  sin(Oy)  in  E,  where  8  is  the 
approximate  "frequency"  of  u  in  E,  then 


and  by  (5.6),  since  C(p)  =  (p!) 


-2 


for  Taylor's  theorem, 


(5.13) 


tE  =  c2  (pi)-2  h2  e2(p+1)  h2p 
=  c2  (e  h)2(p+1>  (pi)-2. 

Since  0  is  independent  of  h  end  p,  by  (5.13), 

o  -  (p  tE  1  h)V1  h_1-  (5.14) 

Using  the  erode  bound  d(pl)/dp  £  p  pl/2, 

3(5E)/9p  =  ?E  (2  In (0  h)  -  p), 

so  we  estimate 

' 

4P>  “  <  ?E  (“2  ln(«  h)  ♦  p)  /  w»  if  p  <  m-1  (5.15) 

if  P  l  m-1. 

In  practice,  we  estimate  by  the  squared  error  indicator  1 1  tall  I 
(or  I  MalM2  ,  MUMI2  ,  or  s2)  in  formulae  (5.12),  (5.14),  (5.15),  and 

AS 

the  formulae  for  pg  . 

Expressions  (5.12)  and  (5.15)  indicate  that  hg  tends  to  decrease  as 
E  is  refined.  Thus,  the  adaptive  refinement  algorithm  tends  to 
approximately  equalise  {Xg}.  When  the  order  is  uniform,  this  reduces  to 
"asymptotic  equidistribution"  of  the  errors  among  the  elements  (e.g.. 


Fried  and  Tang  [20]) 


w  N 
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5.4  A  symmetrised  model  of  predicted  error  behavior 

In  this  section,  we  construct  a  second  model  of  predicted  error 
behavior  based  on  symmetrized  a  priori  error  estimates. 

In  considering  combined  h  and  p  refinement,  Babufka  and  Dorr  [2] 
present  an  a  priori  error  estimate  of  the  form 

1 1 1  el  1 1 2  i  C(kmax)  /£a(z)  G(z,m(z))  h(z)2^(z)p(z)_2<”(z)-1)dz.(5.16) 

Since  C  and  G  are  independent  of  p,  this  estimate  is  more  symmetric  in  h 
and  p  than  (5.5).  Moreover,  it  indicates  the  manner  in  which 
convergence  occurs  when  p  >  m-1 .  Theorems  and  numerical  evidence  in 
[11]  show  the  exponent  in  p  is  optimal  in  the  presence  of  singularities 
of  u  not  located  at  mesh  nodes,  and  that  if  all  the  singularities  of  u 
are  located  at  mesh  nodes,  the  exponent  of  p  can  double.  The  effect  of 
this  doubling  will  be  considered  at  the  end  of  this  section. 

Te  make  the  heuristic  assumption  that  (5.16)  holds  in  each  element, 
with  equality; 

1 1 1 e II 1 2  =  (g  (5.17) 

:*  C(kmax)  /^a(z)  G(z,m(z))  h(x)2p*Z*  p(z)  2^"*Z*  **  dx. 

Assumptions  (5.4),  (5.7),  (5.8),  and  (5.17)  comprise  our  second 
heuristic  model  of  predicted  error  behavior.  The  resulting  estimates  of 
(Xg*^J,  (XgP^ } ,  and  (Xg)  can  be  used  in  the  refinement  algorithm  in 


Seotion  5.3 
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In  our  second  nodcl ,  if  the  order  in  unrefined  element  B  has  never 
been  refined,  then  by  (S.17)  end  (5.7),  we  estimate 

p  ^  (  ln(«f(E))  ~  ln<*E>  )  /  ln<4)  -  1. 

On  the  other  hand,  if  the  order  in  B  has  been  refined,  then  if  p  <  m-1, 

?B  /  ?"  -  h2  (p/(p-i)>‘2<*"1), 

else  if  p  2  »-l# 

?E  /  tE  “  (P/(P-1))"2<1,_1). 

Since  the  function  (p/(p-l))  2*p-1^  equals  1/4  at  2,  and  is  monotone 

-2 

decreasing  with  limit  e  as  k  goes  to  «*,  ve  decide  whether  to  increase 

M 

our  estimate  of  p  on  the  basis  of  (5.11),  vith  pCDT  ■  a  .In  order  for 
our  estimate  of  p  to  be  logically  consistent,  ve  restrict  it  to  the 
interval  (0,pg). 

As  in  Section  5.3,  ve  estimate 

k*,h)  -  (  p  /  w  )  .  (5.18) 

If  p  <  m-1,  ve  estimate 

4P)  -  tB  <“*  i»(h)  +  2  (m-l)/p)  /  v' .  (5.18) 

Otherwise,  ve  estimate 


4P)  -  <E  (m-l)/p)  /  v' 


(5.20) 
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In  practice,  we  approximate  5g  by  the  squared  error  indicator 
lllelllg  (or  IMelMg  ,  IHsIMg  ,  or  8g)  in  formulae  (S.18),  (5.19), 
(5.20),  and  the  formulae  for  p. 

Suppose  the  exponent  of  p  in  (5.17)  is  doubled.  Then  our  estimate 
of  Xg^  is  exactly  as  before,  and  our  estimate  of  XgP^  changes  by  a 
factor  of  2  when  p  2  m~l,  and  by  a  small  factor  when  p  <  m-1  and  hg  is 
small.  The  main  effect  is  to  encourage  order  refinement  when  m  is  small 
and  p  >  m-1.  But  then  Xgp^  is  usually  much  smaller  than  X^^ ,  so  the 
practical  effect  on  refinement  behavior  is  slight.  In  general,  the  form 
of  (5.17)  is  more  important  than  the  values  of  the  exponents. 

5.5  Choosing  the  refinement  omtoff 
A  simple  choice  for  X^^.  is 

XCUT  “  C 

for  some  constant  0  <  C  <  1,  where 

X-  -  max{A_) . 

*  E  E 

If  it  too  large,  (5.2)  may  not  hold.  If  X^^  is  too  small,  we 

may  refine  an  element  with  small  Xg  when  it  may  be  more  efficient  to 
refine  an  element  with  large  Xg  twice.  Ve  can  avoid  this  situation  by 
letting  Xgjpj,  approximate  the  largest  Xg  for  a  once-refined  element. 

This  quantity  can  be  estimated  by  the  marginal  efficiency  of  refinement 


for  element  B  after  it  has  been  refined  once: 
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If 


and  the  order  in  B  has  never  been  refined#  we  can  estivate 


~  ~  2 
XCUT  (  pg  /  WE  )  1  ?f(B)  ‘ 

If  Xjj  “  and  the  order  ia  B  has  been  refined,  then  by  (5.13)  and 

(5.12)  (or  (5.17)  in  onr  second  model),  we  can  estimate 


l  =  4"(pB+1)  klh) 

W  4  E  kE 


(P > 


If  Kg  =  .j  ,  then  in  onr  first  model,  by  (5.13)  and  (5.15),  we  can 
estimate 


xctrr  :  4"’  <Wl  '  'V1’  <  *,<V,,'W'(V1>  >• 

and  in  onr  second  model,  by  (5.17)-(5 .20) ,  we  can  estimate 
*COT  =  XEP>  4  <<PE+1>'^1>~2<",I"1,  <  *'(PB)/w'(PB+l)  ). 


5.6  Drawbacks 

The  heuristic  assumptions  leading  to  onr  estimates  of  4^  a*y 

too  crude  to  yield  useful  refinement  information.  We  base  an  alternate 

( 

refinement  algorithm  solely  on  Xg  and  p: 
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1.  Find  the  unrefined  eleaent  B  with  the  largest  X^ . 

£ 

2.  Determine  X^  <  Xgh) . 

3.  For  each  unrefined  element  E  with  xi1^  2  X, .  : 

£  cut 

a.  If  p  >  p,  refine  the  mesh  in  E. 

b.  If  p  i  p,  refine  the  order  in  E. 


4.  Take  to  be  the  resulting  finite  eleaent  space. 

If  the  order  in  element  E  is  refined,  p  is  estimated  by  (5.11). 
Then  if  and  £g  behave  similarly  to  I  I le(p-l) 1 1  lg  and  I  I  lei  I  lg  ,  p 
tends  to  be  one  more  than  its  (usually)  optimal  value  m-1.  Conversely, 
if  tg  <<  1 1 1 e(p-l) 1 1 1 g  ,  then  p  tends  to  be  less  than  m-1.  The  latter 
situation  arises  most  frequently  when  kg  >  kg  on  the  sides  of  E.  Thus, 
in  our  current  code,  if  kg  2  kg+1  on  any  three  sides  of  E,  or  kg  2  kg+2 
on  any  side  of  E,  then  the  next  time  E  is  refined,  the  order  is  refined, 
regardless  of  the  size  of  p. 

Lack  of  de-refinement  penalizes  early  order  refinement.  For 
example,  suppose  h  i  1/4  and  k  2  3  are  desired  near  the  origin.  If  the 
order  is  refined  first,  then  k  2  3  everywhere  (Figure  5-1). 

The  resulting  finite  eleaent  spaoe  may  not  satisfy  (5.3),  even  if 
the  heuristic  models  of  error  behavior  hold.  For  instance,  if 
X^^  <  XgP) ,  it  may  still  be  more  efficient  to  refine  B  in  h  if  the 
resulting  decrease  in  {g  is  greater  with  mesh  refinement  than  with  order 
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Figure  3-1:  Early  order  refiaeaeat  peaalty 
refiaeaeat.  sad  if  the  other  hg's  are  very  saall. 

Rigorous  aodels  are  preferable  to  the  heuristic  aodels  used  ia  this 
ohspter.  A  rigorous  approaoh  for  oae-diaeasioaal  problems  is  developed 
ia  Babulka  sad  Rheiaboldt  [7]. 

3.7  Expected  error  behavior 

Ia  this  sectiou,  we  aodel  the  asyaptotie  behavior  of  error  with 
respect  to  work  for  three  problea  classes  consisting  of  probleas  with: 

1.  saooth  solutions. 

2.  solutions  with  point  singularities. 

3.  solutions  with  line  singularities. 
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We  use  the  heuristic  model  from  Section  5.3.  We  expect  similar  behavior 
using  the  heuristic  model  from  Section  5.4. 

5.7.1  A  smooth  solution 

Suppose  u  is  smooth  everywhere  in  Q,  and  h(z)  and  p(z)  are  uniform. 
By  assumption  (5.6), 

S  =  C(p)h2p, 

and  by  assumption  (5.4), 

W  =  pP  h”2.  p  2  2. 

Minimizing  £  while  holding  W  fixed,  2  p/h  £  -  X.  2/h  W  =  0  and 

((C'/C)  +  2  ln(h))£  +  A  0/p  W  ®  0  imply 

t  ~  -r/2 
h  =  e  ' 

and 

?  =  C  e'rP. 

where 

y  :=  P  +  C'/C. 

Thus,  it  is  best  to  fix  the  mesh  and  only  refine  the  order.  Suppose 
C(p)  -  AP  ,  A  >  1. 

Since  C'/C  =  ln(A) ,  then  A  e  *  =  A  e  ^  <  l,  so 

£  =  (A  e’y)p 

is  exponentially  decreasing  in  p.  But  W  *  p^  e^  is  polynomial  in  p,  so 
if  there  exists  A  <  •  such  that  C(p)  £  Ap  for  sufficiently  large  p,  then 
£  is  exponentially  decreasing  with  respect  to  W. 
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5.7.2  A  point  singularity 

Suppose  u  has  an  isolated  singularity  at  a  point  a  in  Q,  such  that 

2 

near  a,  u  has  1+a  derivatives  in  L  ,  and  u  is  saiooth  away  from  or. 
Suppose  there  are  unrefined  elements  at  levels  1  to  L  in  the  mesh  for 
S^,  and  at  each  level  only  the  elements  which  touch  the  point  a  are 
originally  refined.  Suppose  p(z)  uniformly  equals  L.  Then 

=  C(a)  h^+2a  if  E  touches  o 

and 

=  C(p)  h^+2p  if  E  does  not  touch  o. 

By  the  proof  of  (2.4),  there  are  no  more  than  a  constant  number  of 
elements  at  each  level,  so 

*  “  5  E  *E 

i  C^p)  {(2_L)2+2a  +  }  <2_i)2+P> 

i  Cj(L)  2_L, 

and 


If  Cj(L)  £  AL  for  some  A  <  2  as  L  goes  to  *,  then  $  is  exponentially 
decreasing  with  respect  to  V.  Babulks  snd  Dorr  [2]  make  a  similar 
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conjecture,  and  under  suitable  assumptions,  prove  that  convergence  is 
faster  than  polynomial. 

5.7.3  A  line  singularity 

Suppose  u  has  a  line  singularity  at  a  line  o  in  Q,  such  that  near 

2 

o,  u  has  1+a  derivatives  in  L  ,  and  u  is  smooth  away  from  or.  To  fix 

ideas,  suppose  p(z)  is  uniform,  a  traverses  Q  from  top  to  bottom,  and 

there  are  elements  at  levels  1  to  L  in  the  mesh  for  S  ,  where  at  each 

level  the  elements  which  touch  the  line  o  are  present  in  S.  .  Then 

L 

=  C(a)  h^+2°  if  E  touches  a 

and 

?E  * 

There  are  at  least  2^  elements  with  level  L  touching  a,  so 

5  l  CjU)  2L  (2~L)2+2a 

and 

W  l  C2<a)  2L. 

Then 

{  1  C3(u)  W"1_2a 

regardless  of  how  p  is  chosen.  Convergence  cannot  be  faster  than 
polynomial,  since  the  complexity  of  geometrically  approximating  a 
asymptotically  dominates  the  cost. 


C(p) 


h2+2P 


if  E  does  not  touch  o. 
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Computational  aspects 


6.1  Introduction 

In  this  chapter,  ve  discus*  some  of  the  computational  aspects  of 
the  methods  presented  in  Chapters  2-5.  Since  ve  have  not  yet 
implemented  all  the  procedures  discussed,  detailed  operation  counts  are 
omitted. 

In  Section  6.2.  ve  consider  the  efficient  assembly  of  the  finite 
element  system  (cf.  Eisenstat  and  Schnltz  [18].  Veiser,  Eisenatat,  and 
Schultz  [31]).  In  Section  6.3,  ve  investigate  the  complexity  of  sparse 
direct  elimination  vith  nested-dissection-type  ordering  of  the  unknovns 
for  the  three  classes  of  problems  discussed  in  Section  5.7.  The 
operation  counts  and  storage  for  the  problems  vith  singularities  are 
linear  (optimal-order)  in  the  number  of  elements.  In  Section  6.4,  ve 
consider  the  efficient  computation  of  the  a  posteriori  error  estimators 

Ve  nov  outline  the  computational  tasks  to  be  performed. 

There  are  four  computational  steps  to  perform  for  eaoh  finite 
element  space  S  in  the  sequenee  (S^) . 
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1.  Fora  the  linear  eye tea  AD*  f,  where 
{B^}  ie  the  beeic  for  S, 

U(x,y)  =  ^  t  DTBT(r,y), 

AIJ  *E(BI,BJ*  ' 

(f.B,)E  . 

After  initializing  A  and  f  to  zero,  for  each  E, 

~E  ~E 

a.  Aeeeable  A  and  f  ,  the  eleaent  etiffneee  aatriz  and 
right  hand  aide  for  the  baeie  (B^),  where  {B^}  ie  the 
baeie  which  would  result  if  E  had  no  irregular  corners. 

E  E  BNjg 

b.  Obtain  A  and  f  froa  A  and  i  ,  by  changing  variables 
froa  {Bj}  to  {Bg} . 

c.  Add  AE  and  fB  to  A  and  f. 


“  2  E  AIJ  ' 

fi  -5Ef?  • 


2.  Solve  A  U  ■  f, 


3.  For  each  B,  coapnte  lllelllg  ,  MleMlg  ,  MleMlg  ,  or  . 

4.  Depending  on  the  error  indicator e  eoapnted  in  Step  3,  either 
stop,  or  adaptively  choose  S^  ,  and  go  to  Step  1. 


Each  integral  v(x,y)  dx  dy  is  approziaated  by 

T  0(E)  T  0(E)  ,  . 

2  s-1  2  t-1  Ws  Wt  v(*,'yt>' 


(6.1) 


whera  and  (wt*yt>  are  the  weights  and  poiats  for  the  G(E) -point 

V 

respectively.  Choices  for  G(E)  are  specified  in  Chapter  7. 

6.2  Assembling  the  linear  system 

Consider  Step  l.a  (Section  6).  For  simplicity*  s oppose 

aE(v,w)  =  (a(r.y)  v.w)E  , 

the  bilinear  form  for  a  weighted  least  squares  problem*  since 
considerations  for  general  are  almost  identical. 

Since  [Bjlg]  is  a  subset  of 

(p®(z)pB(y)  :  i-l*...*kmsz  ,  m-1 * . . . * kmaz  ), 

1  B 

sparse  versions  of  the  tensor-produo t  assembly  procedures  in  Eisenstat 
and  Schultz  [18]  and  Teiser*  Eisenstat*  and  Schultz  [31]  are 
appropriate.  If 

®I  “  pi(x)pm<y>' 

Bj  -  Pj(z)pn(y), 

then*  depending  on  whether  a(z*y)  and  f(z»y)  are  non-separable* 
separable,  or  constant*  A  and  f  can  be  computed  using  the  formulae 


V  ,nd 


Gauss-Legendre  integration  rules  on  [zg. 
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<.<x,y>  Bi.Bj)e  I 

1  t  <WtP«(yt)pn(yt)>  1  1  ,  <Vi(,.)pj(li)>  •(x,*yt) 

(*x(x)«y(y)  Bj.BjJg  - 

1 }  .  <Vitx.)pj(x,)>  S(X.)1  [ }  t  <V.(yt)p«(yt)>  V(yt)1' 


(5I’Ve  "  <  (L  Vi(x.)pj(x.))  <5  t  V«(yt,pn<yt>>  >* 

(fCx.yJ.Bjig  -  J  t  <*tPB<yt>>  I  J  t  <w#p1(xb)>  1* 

(£x(x)fy(y).BI)E- 

[  }  .  <Vi<X,)>  fx(x.)  1  C  2  t  <V«(yt)>  fy(yt)  ]' 
(1'5i)e  =  <  (5  .  Vi<X.,)  (5  t  V-(yt)>  >* 


Hi*  quantities  in  angle  brackets  <>  can  be  precomputed  independent  of  B. 
The  quantities  in  square  braokets  tl  should  be  eoaputed  as  intermediate 
suae  in  each  E.  Since 


A8 

AIJ 


HS 

A(i,m)(j,n) 


A( i,n) ( j ,m) 


A8 

( j»n) ( i»n) 


A(J.n)(i,m)  * 
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i,m)(j,n)  " 

T  rE 

L  L  VL  I  i,L) ( j ,n) 


.E 


(j,n)(i»n) 

\  P  Vn.P  [  \  L  %.L  X(i.L)<j,P> 


and 


P?.  ,  -  <  f?,  *  if2. 

(i.m)  1  (i.m) 

\  L  Vn.L  *<i.L)  i-2* 


if  2,  jf2, 
i-2.  Jf2, 
if 2,  j-2. 


i-2,  j-2. 


The  resulting  computation  is  relatively  expensive  when  kg  is  small. 

When  a(x,y)  is  separable,  fewer  mnltipliea  are  required  if  Step  l.b  is 
E  E 

omitted,  and  A  and  f  are  formed  directly. 


(.3  Solving  the  linear  system 

Consider  Step  2  (Section  6),  solving  the  linear  system  A  0  -  f.  In 
this  section,  we  investigate  the  complexity  of  sparse  direct  elimination 
with  nested-dissection-type  ordering  of  the  unknowns  for  the  three 
classes  of  problems  discussed  in  discussed  in  Section  5.7. 

<•3.1  A  smooth  solution 

As  in  Section  5.7,  suppose  U  is  smooth  everywhere  in  0.  Suppose  k 

-2 

and  h  are  uniform,  and  there  are  N  :•  h  elements.  There  are 
(k-4)(k-3)/2  element  basis  functions  nonxero  in  each  E.  In  sparse 


Gaussian  elimination,  if  the  element  baaia  fnnctiona  ere  ordered  firat, 
their  elimination  indnoea  no  fill-in  (i.e.,  they  are  ataticellw 
condenaed  [28]).  This  initial  elimination  requires  -  N((k-4)(k-3)/2)*/6 
multiplies  and  *  N((k-4) (k-3)/2)^/2  storage.  Although  this  coat 
dominates  when  k  >>  N,  in  solving  problems  to  usual  aeouraoy, 
elimination  of  the  element  basis  functions  is  relatively  inexpensive. 
There  are  -  (2k-3)N  remaining  basis  functions  (k-2  side  basis  functions 
associated  with  each  distinct  element  side,  and  1  vertex  basis  function 
associated  with  each  distinot  vertex). 

Consider  a  nested  disseetion  ordering  of  the  side  and  vertex  basis 
functions,  in  conjunction  with  sparse  elimination  of  A  (e.g.,  George  and 
Liu  [22]).  A  sample  nested  dissection  ordering  is  shown  in  Figure  <-2. 


Figure  €-2 j  Nested  disseetion  ordering  with  a  smooth  solution 
The  unknowns  are  ordered  according  to  the  groups  numbered  l-d.  Groups 
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1-4  correspond  to  basis  functions  on  uniform  snbmeshes  with  N/4 
elements:  in  each  group  the  unknowns  are  recursively  ordered  using 
nested  dissection  on  the  subnesh.  Since  group  5  separates  the  unknowns 
in  the  firs^t  four  groups,  the  number  of  multiplies  to  eliminate  the 
unknowns  in  groups  1-5  is,  to  dominant-order  terms,  bounded  by  M(N) , 
where 

M(N)  -  4  M(N/4)  +  C  (IN1'2)3. 

The  rightmost  term  is  a  bound  on  the  number  of  multiplies  needed  to 
densely  eliminate  the  unknowns  in  group  5.  Hie  solution  to  this 
recurrence  relation  is  0(kV/2>  (e.g..  Rose  and  Whitten  [26],  Theorem 
1),  as  is  the  operation  count  to  eliminate  the  remaining  unknowns  in 
group  6.  Thus,  the  expected  operation  count,  neglecting  static 
condensation  cost,  is  0(kV/2).  Similarly,  the  expected  storage  is 
0(k2Nlog(k2N)). 

6.3.2  A  point  singularity 

As  in  Section  5.7,  suppose  u  has  an  isolated  singularity  at  the 
point  o*(0,0) .  Suppose  there  are  unrefined  elements  at  levels  1  to  L  in 
the  mesh,  and  at  each  level  only  the  elements  which  touch  e  are  present. 
Suppose  k  is  uniform. 

There  are  N  ■  3L  unrefined  elements  in  the  mesh.  After  statie 
condensation,  “  (2k-3)N  side  and  vertex  basis  functions  remain. 

Consider  a  recursive  ordering  of  the  side  and  vertex  basis 
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functions,  in  conjunction  with  sparse  elimination  of  A.  A  saaple 
ordering  is  shown  in  Fignre  6-3. 


Figure  6-3:  Recursive  ordering  with  a  point  singularity 
The  unknowns  are  ordered  aooording  to  the  groups  nusibered  1-3.  Group  1 
corresponds  to  basis  functions  on  a  subnosh  with  N-3  elements:  the 
unknowns  in  group  1  are  recursively  ordered  on  the  submesh.  The  number 
of  multiplies  to  eliminate  the  unknowns  in  group*  1  and  2  is,  to 
dominant-order  terms,  bounded  by  M(N) ,  where 

M(N)  -  M(N-3)  +  C  k3. 

The  rightmost  term  is  a  bound  on  the  number  of  multiplies  needed  to 

densely  eliminate  the  unknown*  in  group  2.  The  solution  to  this 

a 

recurrenoe  relation  is  0(k  N) ,  and  the  operation  count  to  eliminate  the 
remaining  unknowns  in  group  3  is  O(k^).  Thus,  the  expected  operation 
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a 

const,  neglecting  static  condensation  cost,  is  0(k  M) .  Similarly,  the 

2 

expected  storage  is  0(k  N) . 

(.3.3  A  line  singularity 

As  in  Section  5.7,  suppose  u  has  a  line  singularity  along  the  line 
a  ■  { (O.y) ) .  Suppose  there  are  unrefined  elenents  at  levels  1  to  L  in 
the  aesh,  where  at  each  level  the  elements  which  touch  the  line  a  are 
present.  Suppose  k  is  uniform. 

M  L 

There  are  N  -  3(2  )  unrefined  elements  in  the  mesh.  After  static 
condensation,  =  (2k-3)N  side  and  vertex  basis  functions  remain. 

Consider  a  nested-dissection-type  ordering  of  the  side  and  vertex 
basis  functions,  in  conjunction  with  sparse  elimination  of  A.  A  sample 
ordering  is  shown  in  Figure  (-4.  The  unknowns  are  ordered  according  to 
the  groups  numbered  1-4.  Groups  1  and  2  correspond  to  basis  functions 
on  submeshes  with  (N-2)/2  elements:  in  each  group  the  unknowns  are 
recursively  ordered  on  the  subaesh.  Since  group  3  separates  the 
unknowns  in  the  first  two  groups,  the  number  of  multiplies  to  eliminate 
the  unknowns  in  groups  1-3  is,  to  dominant-order  terms,  bounded  by  N(N) , 
where 

M(N)  -  2  M(N/2)  +  C  (klogN)3. 

The  rightmost  term  is  a  bound  on  the  number  of  multiplies  needed  to 

densely  eliminate  the  unknowns  in  group  3.  The  solution  to  this 

3 

recurrence  relation  is  0(k  N)  (e.g..  Rose  and  Whitten  [2(]),  and  the 
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Figure  6-4:  Nested-dissection-type  ordering  with  e  line  singularity 

operation  count  to  eliminate  the  remaining  unknowns  in  group  4  is 
3 

0( (klogN)  ).  Thus,  the  expected  operation  count,  neglecting  static 

3 

condensation  cost,  is  0(k  N) .  Similarly,  the  expected  storage  is 
0(k2N) . 

6.3.4  Discussion 

The  operation  counts  and  storage  for  the  problems  with 
singularities  are  optimal-order  in  the  number  of  elements  N.  Ve  expect 
similar  operation  counts  for  more  general  problems  with  point  or  line 
singularities.  The  optimality  for  the  line  singularity  problem  depends 
on  the  fact  that  we  can  systematically  find  sets  of  oo^3-*)  elements 
(e  >  0)  which  separate  the  remaining  elements  into  two  or  more  disjoint 
sets  of  equal  size. 
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We  do  not  know  an  effioient  (i.e.,  0(N)),  automatic,  way  of  finding 
oxderings  such  as  those  presented.  See  George  and  Liu  [22]  ,  Sections 
7.2  and  8.2,  for  work  in  this  direction. 

In  the  example  problems,  the  orderings  of  the  interior  unknowns 
correspond  to  the  reverse  of  the  order  in  which  the  edges  associated 
with  the  unknowns  are  created.  This  correspondence  does  not  always 
hold.  (For  example,  suppose  u  has  a  point  singularity  at  (1/2,0).) 

6.4  Computing  the  error  indicators 

Consider  Step  3  (Section  6),  the  computation  of  the  error 
indicators  lllelllg  ,  IjJsMIg  *  IJJslilg  '  or  =E  * 

F£(v)  (f.v)E  +  <[a3U/dn]  ,v>aE  -  SgOJ.v). 

First,  suppose  E  has  no  irregular  corners.  Then  the  possible  quantities 
to  be  computed  in  Step  3  are 

1.  (aE(BI,BJ)} ,  where  {Bj}  is  the  basis  for  Sg  which  would 
result  if  E  had  no  irregular  corners.  These  quantities  can 
be  computed  as  in  Step  l.a. 

2.  {(f,B|.)E  :  §j  in  Sg-Sg) .  These  quantities  can  be  computed 
as  in  Step  l.a. 

3.  (agdJ.Bj)  :  Bj  in  Sg-Sg).  First,  convert  the  representation 

0(x.r)lE- J  ,  o,  “,(,.7) 
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into  the  form 

For  example,  for  element  E  in  Figure  6-1,  when  k  is  uniform, 
if 


'’“•’"'e-LwLv.) 

2  m  0(2,m) 


Pj(x)  p®(y)  + 

1  ID 

E,  .  f (E) .  . 

p 2<*>  pb  (y>* 


then 


t  i,m) 


=  /u 


( i,m) 

\  L»L(m)  VL,m  ^(i.L) 


where 


i  *  2, 

i  *  2, 


L(m)  = 


m  #  2, 

m  =  2, 


Other  cases  are  handled  analogously. 

Then,  for  lllelllg  ,  compute  (a^(U,B^)}  with  a  matrix-vector 
multiply.  For  NleMlg  or  IMeMlg  ,  compute  (U(xs,yt)}  by 
the  formula 


B<vV  ‘L^t1 1 }  i  pi(x»)  s(i,«) ] 


(6.3) 
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and  then  use  the  formula 
(a(x,y)  U.BjJg  =  J 

1  5  s  >  1  °(*,>yt)*<xg*yt>  1  1 • 

4.  {< [a3U/3n] ,Bj>gE) .  In  a  preprocessing  step,  compute  and 

store  the  coefficients  of  the  one-dimensional  polynomials 

3U/3nlg  on  each  element  side.  Then,  for  example,  on  the  left 

side  of  element  E,  compute  {adU/dn(x*,y.  ) |_}  and 

Etc 

{a3U/8n(xjj,yt) Ig,} ,  where  element  E* ,  if  it  exists,  shares 
part  of  the  left  side  of  E.  Form 

<trtir/i.i  „ft)  : 

5  t  <V.<ytl/2n,atl/i"“<VrtllE  '  ■iWS"<VVlE.) 

for  each  B_  =  p.(x)p  (y) .  Other  cases  are  handled 
1  JL  O 

analogously. 

5.  HleMlg  ,  IMcIMg  »  or  IH*Mlg  •  Suppose  the  system  for 

The  storage  required  can  be  reduced  by  intermixing  the  preprocessing 
step  with  error  indicator  computation. 
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7.1  Introduction 

In  this  chapter,  ve  present  nnsierical  results  obtained  using 
prototype  codes  written  in  FOSHAN,  which  iaplenent  some  of  the  nethods 
discussed  in  Chapters  2  -  S.  The  problea  set,  containing  one  problen 
from  each  of  the  three  classes  considered  in  Section  5.7,  is  specified 
in  Section  7.2.  In  Section  7.3,  we  present  the  error  behavior,  which  is 
in  fair  accord  with  Section  5.7.  In  Section  7.4,  we  present  the  error 
estiaator  behavior:  the  error  estiastors  are  usually  accurate  to  within 
a  factor  of  two.  In  soae  oases,  soae  of  the  estiastors  appear  to 
converge  to  the  nora  of  the  true  error.  In  Section  7.5,  we  present  aore 
nnaerical  evidence  of  error  estiaator  convergence. 

7.2  The  problea  set 

Prototype  codes  were  written  in  FORTRAN,  iaplenent ing  soae  of  the 
aethods  disoussed  in  Chapters  2-5.  Nuaerioal  tests  were  conducted  to 
investigate  the  behavior  of  the  aethods. 

The  test  problea  set  consisted  of  three  Poisson  probleas  with 


Dirichlet  boudtry  conditions 


2  2 

-  D  u  -  D  u  ■  f(x,y)  in  the  interior  of  Q, 

*  y 

n(x,y)  *  g(x,y)  on  dQ, 


with  eolations  specified  in  Table  7-1. 


Problen  I 

1  Solution  class 

1 1  Solution 
_i  i 

1  1 

1  snooth 

1  1 

11 

1 1  u(x,y)  ■  cos(e  +xy) 

II 
.1 1 

2  1 

1  point 

II 

II 

1 1  u(x.y)  * 

1 1 

r1/2  sin(9/2)  , 

i  singularity 

1  • 

1 1  x  ■  r  cos (6)  ,  y  =  r  sin (6) 

II 
.1 1 

3  ! 

1  line 

1  1 

1 1  u(x,y)  -  i 

n 

H 

1 

N't 

+» 

O 

K _ 

1  singularity 

jj 

1 1 

sin(nt)  -1/2  <  t  <  1/2 

1 1 
ii 

M 

Jl  1/2  £  t  . 

s 

■ 

lit*  5x-y-l 

s  • 

Table  7-1:  Test  problems 

The  solution  to  Problen  1  is  in  H^O)  for  ell  k.  The  solution  to 
Problen  2  is  the  fundenentsl  eigenfunction  for  s  crsok  problen  with 
interior  angle  2n;  this  function  is  in  H 3/2“«  '(0)  for  any  s  >  0,  but  not 
in  The  solution  to  Problen  3  is  a  sneered  shock  ware  with 

shock-width  1/5  in  x;  this  function  is  in  but  not  in  for 
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any  a  >  0. 

Ve  tasted  the  seven  refinement  methods  specified  in  Table  7-2.  All 
adaptive  sequences  of  finite  element  spaces  started  out  with  2X2  meahes. 
The  adaptive  sequences  for  methods  +  and  X  atarted  out  with  k_  ■  2. 


II 

1 1 

Method 

!  4 

II 

‘e 

II 

1  1 
II 

1 1 

0 

1  uniform  (1/2, ... ,1/32) 

1  . 

II 

-  2 

II 

1  1 
II 

1 1 

X 

1  -  1/2 

1 

II 

adaptive  (Section  5.4) 

II 

1  1 
II 

1 1 

2 

1  adaptive  (Section  5.3) 

if 

-  2 

II 

1  1 
II 

1 1 

3 

1  adaptive  (Section  5.3) 

II 

-  3 

II 

1  1 
II 

4 

1  '  '  '  '  1  1  r 

1  adaptive  (Section  5.3) 

’ll" 

-  4 

II 

II 

+ 

1  . 

1  adaptive 

’ll’ 

adaptive  (Section  5.6) 

II 

II 

II 

X 

I  adaptive 

’ll’ 

II 

adaptive  (Section  5.4) 

II 

II 

Table  7-2*  Refinement  methods 

We  used  three  heuristic  refinement  procedures,  based  on  the 
formulae  in  Chapter  5.  In  the  "true  error"  procedure,  ve  aet 
{g  ■  I  Mel llE  .  In  the  "approximate  error"  procedure,  ve  set 
tg  ”  I  Ms  I  Mg  ,  vith  S  composed  of  piecevise  polynomials  one  order 

2 

higher  than  S.  In  the  "residuals  and  jumps"  procedure,  ve  set  (g  ■  Sg 

EST 

For  simplicity,  ve  set  Xgpf  -  To  avoid  ill-oonditioned  linear 

systems,  ve  restricted  kg  £  7.  Our  tests  vers  run  in  single  precision 
on  a  DEC-Systen  2060  computer. 
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G(E)-point  Gauss-Legendre  quadrature  rulea  (6.1)  were  used  to 
approximate  integrals  in  each  element  E.  In  the  original  assembly 
phase,  G(E)  =  kmaxg  ,  where 

kmaxg  :«  max {kg,  :  E'  is  in  N(E>). 

In  the  error  indioator  phaae,  G(E)  »  kaaz^+1.  In  the  evaluation  of 
lllelllg,  G(E)  “  kmaXg+3 .  With  these  values,  for  problems  with  smooth 
solutions,  the  errors  due  to  numerical  quadrature  are  asymptotically 
negligible  (e.g.,  Ciarlet  [17],  Section  4.1). 

7.3  Error  behavior 

In  this  section,  the  resulting  errors  for  the  adaptive  aequencea  of 
finite  element  spaces  are  depicted  on  lo*10-1°*10  •cale,‘  Th® 
x-coordinate  of  each  data  point  is  the  number  of  nonzeroes  in  the 
stiffness  matrix  ,  and  the  jr-coordinate  is  I  llelll.  Since  the  true  work 
generally  grows  faster  in  k  than  the  number  of  nonzeroes,  the  measure  of 
work  is  biased  in  favor  of  high-order  methods. 

Figure  7-1  presents  the  errors  due  to  applying  the  "true  error" 


e 

Our  prototype  codes  should  be  made  more  efficient  before  computer 
time  is  used  as  the  measure  of  work. 


procedure  to  Problem  1.  Sinoe  u(x,y)  la  smooth,  mulleta  uik 

refinement  does  not  improve  the  asymptotic  convergence  rate.  The  alopea 

Jr_l  4 

of  tbe  lines  correspond  to  0(h  )  error  and  0(h  k  )  vork.  Since  k 

increases  for  methods  E,  +,  and  X,  these  methods  aahieve  faster  than 
polynomial  convergence.  Methods  +  and  X  perform  some  initial  mesh 
rfinement  (see  Figure  7-2,  with  kg  labelled  in  each  element),  so  they 
are  somewhat  less  efficient  than  method  K. 

Figures  7-3  and  7-4  present  the  errors  dne  to  applying  the 
"approximate  error"  and  "residuals  and  jumps"  .cednxos  to  Problem  1. 
The  errors  are  similar  to  the  ones  shown  in  Fig^:  »  7-1, 

Figure  7-5  presents  the  errors  due  to  applying  the  "true  error" 

procedure  to  Problem  2.  The  slope  of  the  line  for  method  V  corresponds 
1/2  —2 

to  an  0(h  )  convergence  rate  and  an  0(h  )  work  count.  Method  K  has 

almost  the  same  line,  with  slope  corresponding  to  an  0(k  *)  oonvergenee 
4  e 

rate  and  an  0(k  )  work  oount.  The  slopes  of  the  lines  for  methods  2,  3, 

e 

Vith  the  work  oount  proportional  to  the  number  of  unknowns,  method  E 
would  appear  twice  as  efficient  as  method  D  (see  [11]).  Vith  a  more 
realistic  work  count,  method  E  would  asymptotioslly  be  the  least 
efficient  method  considered. 
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and  4  are  the  same  as  the  slopes  for  the  sane  Methods  in  Pignre  7-1. 
Local  nesh  refinement,  in  effeet,  transforms  variables  so  the  problem 
looks  like  it  has  a  smooth  solution  on  a  uniform  mesh.  Methods  +  and  X 
are  fairly  efficient  in  all  cases,  refining  the  mesh  near  the 
singularity,  and  refining  the  order  where  the  solution  is  smooth  (see 
Figure  7-6) . 

Figures  7-7  and  7-8  present  the  errors  due  to  applying  the 
"approximate  error"  and  "residuals  and  jumps"  procedures  to  Problem  2. 
The  errors  are  similar  to  the  ones  shown  in  Figure  7-5. 

Figure  7-9  presents  the  errors  due  to  applying  the  "true  error" 

procedure  to  Problem  3.  The  slope  of  the  line  for  method  D  corresponds 

to  an  0(h)  convergence  rate  and  an  0(h  )  work  count.  The  line  for 

method  K  appears  to  level  out,  reflecting  an  0(k*)  work  count  and  slow 

convergence  in  k.  The  slope  of  the  line  for  method  2  reflects  an 
-1/2 

improving  (<  N  )  convergence  rate  and  an  0(N)  work  count,  where  N  is 
the  number  of  elements.  At  the  modest  accuracies  required,  the  mesh  is 
refined  not  only  along  the  line  singularities  at  the  edges  of  the 
boundary  layer,  but  in  the  boundary  layer  as  well  (see  Figure  7-10).  Te 
expect  the  slope  of  the  line  for  method  2  to  approach  -1  asymptotically, 
reflecting  an  0(N  *)  convergence  rate.  The  slopes  of  the  lines  for 
methods  3,  4,  +  ,  and  X  refleot  convergence  rates  tending  toward  0(N  1). 
Figure  7-11  depicts  a  sample  grid  for  method  +.  Note  that  the  shape  of 
the  boundary  layer  cannot  be  discerned  from  the  grid. 


107 


Figures  7-12  end  7-13  pretest  the  errors  due  to  applying  the 
"approximate  error"  and  "residuals  and  jumps"  procedures  to  Problem  3. 
The  errors  sre  similar  to  the  ones  shown  in  Figure  7-9. 

In  this  test,  methods  +  and  X  appear  to  be  the  most  robust,  since 
they  are  fairly  efficient  for  all  three  problems,  regardless  of  whether 
the  errors  are  large  or  small.  Eaoh  of  the  other  methods  does  poorly 
for  some  problems  in  some  cases.  Table  7-3  smamarizes  the  observed 
efficiency  a,  where 

1 1  lei  1 1  13  (the  number  of  nonzeroes)  a. 


II 

1 1 

Method 

1  Problem  1  1 

1  Problem  2 

i  . . 

II 

Problem  3 

II 

1  1 
II 

1 1 

D 

1  1/2  | 
j  j 

i  1/4 

if 

< 

1/2 

II 

1  1 
II 

1 1 

K 

1  <k-l)/2  | 

j  j 

1  1/4 

II 

** 

m 

1/2 

’ll 

1 1 
II 

1 

2 

1  1/2  | 
|  l 

1  1/2 

if 

> 

1/2 

II 

1  1 
II 

1 1 

3 

1  1 

1  1  1 

1  1 

1  1 

if 

m 

1 

II 

1  1 
II 

1 

4 

1  3/2  I 

j  j 

1  3/2 

II 

- 

1 

II 

1 1 
II 

1 1 

+ 

1  (k-l)/2  | 

|  [ 

1  (k-l)/2 

II 

** 

1 

II 

1  1 
II 
II 

X 

1  <k-l)/2  I 

1  1 

1  (k-l)/2 

if 

II 

1 

II 

II 

Table  7-3:  Efficiency  summary 
7*4  Error  estimator  behavior 

Figures  7-14  -  7-19  depict  the  relative  errors  p  and  £  in  the  error 
estimators  for  the  tests  run  using  the  "approximate  error"  and 
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"residuals  and  jumps"  procedures  respectively. 

In  the  "approximate  error"  procedure,  the  error  is  usually 
estimated  within  a  factor  of  two  (-.5  ip  <.  1) ,  and  always 

-.7  i  p  i  1.6  . 

p  appears  to  converge  to  zero  in  methods  U  and  2  for  Problem  1,  and  in 
methods  2  and  3  for  Problem  2  (see  Section  7.5).  For  Problem  1,  p 
varies  somewhat  erratically  in  methods  +  and  X,  especially  when  kg  is 
large.  This  may  be  partly  due  to  ill-conditioning. 

In  the  "residuals  and  jumps"  procedure,  the  error  is  always 
estimated  within  a  factor  of  two.  In  particular,  in  method  3, 

-.25  i  g  1  .2°  , 

and  g  appears  to  converge  to  zero.  In  methods  +  and  X,  g  varies  more 
smoothly  than  p.  Thus,  in  this  test,  the  "residuals  and  jumps" 
procedure  appears  to  be  slightly  more  robust  than  the  "approximate 
error"  procedure. 

7.5  Brror  estimator  convergence 

As  further  evidence  of  apparent  error  estimator  convergence,  we  now 
present  the  behavior  of  the  four  error  estimators  illtlll,  1 1.1  a  Ij.  I, 
IMcMI,  and  s,  with  relative  errors  p,  g,  p,  and  g  respectively,  for  a 
sample  variable-coefficient  problem  with  a  smooth  solution,  k  and 
h  “  1/n  are  uniform.  S  oontains  piecewise  polynomials  one  order  higher 
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than  S.  Ve  have  observed  similar  error  estimator  behavior  for  other 
problems  with  smooth  solutions. 

Table  7-4  depicts  the  relative  errors  in  the  error  estimators  for 
the  Dirichlet  problem 

-  D  (ex^D  n)  -  D  (e^D  n)  +  a/(l+x+y)  *  f(x,y) 
xr  y  y 

in  the  interior  of  Q, 

n(x.y)  =  0  on  dQ, 

where  f(x.y)  is  chosen  so  that  a(x,y)  =  exy  sin(nx)  sin(ny) . 

In  Table  7-4,  the  error  estimatora  are  always  within  a  factor  of 
two  of  lllelll.  None  of  the  error  estimators  converge  to  1 1  lei  1 1  with  k 
increasing  and  h  fixed,  bat  some  appear  to  converge  with  h  decreasing 
and  k  fixed.  Illslll  and  I  Mel  M  appear  to  converge  when  k  £  4. 

I  MsIH  appears  to  converge  when  k  <  4,  bat  does  not  converge  when 
k  *  4.  The  convergence  rates  are  not  clear.  When  k  *  2  or  3,  s  behaves 
like  a  constant  times  lllelll. 
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!  ii 

1 1. 

n 

II 
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II 

A* 

P 
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Table  7-4 i  Error  estimator  behavior 


Figure  7-1:  Problem  1  -  true  error 
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Figure  7-3 s  Problem  1  -  approximate  error 
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Figure  7—4:  Problen  1  —  residuals  and  juaips 
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W|*t*  7~5t  Problem  2  -  trw  error 


P0INT  SINGULARITY  -  TRUE  ERR0R 
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Figure  7-7:  Problem  2  -  approximate  error 


P0INT  SINGULARITY  -  APPR0XIMATE  ERR0R 
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Figure  7-8 j  Problea  2  -  residuals  and  jnaps 
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Fifnra  7-11:  Staple  grid  for  problem  3  (method  +.  residual*  and  jumps) 
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Figure  7-12:  Problem  3  -  approximate  error 
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Figure  7-13 1  Problea  3  -  residuals  sad  j naps 
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Figure  7-14:  Problem  1  -  approximate  error 
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Figure  7-15:  Problea  1  -  residual*  and  jnaps 
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Fifttx#  7-16:  Problem  2  —  approximate  error 
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Figure  7-17:  Problem  2  -  residuals  sad  Jaaps 
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Figure  7-18:  Problem  3  -  approximate  error 
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